cont, 2

docs/src/man/expectedCFs.md

@@ -35,6 +35,7 @@ R"plot(0:1, 0:1, type='l', bty='L', lwd=0.3, col='#008080', xlab='quartet CF obs
 R"set.seed"(1234); # hide
 R"points(jitter($obsCF,amount=0.005),jitter($expCF,amount=0.005),col='#008080',bg='#00808090',pch=21)"; # hide
 R"dev.off()"; # hide
+nothing # hide
 ```
 To install ggplot2 if not installed already, do:
 `R"install.packages('ggplot2', dep=TRUE)"`

docs/src/man/snaq_plot.md

@@ -47,6 +47,7 @@ R"svg(name('snaqplot_net0_1.svg'), width=4, height=3)" # hide
 R"par"(mar=[0,0,0,0]) # hide
 plot(net0, :R);
 R"dev.off()"; # hide
+nothing # hide
 ```
 ![net0_1](../assets/figures/snaqplot_net0_1.svg)
 
@@ -69,6 +70,7 @@ R"svg(name('snaqplot_net1_1.svg'), width=4, height=3)" # hide
 R"par"(mar=[0,0,0,0]) # hide
 plot(net1, :R, showGamma=true);
 R"dev.off()"; # hide
+nothing # hide
 ```
 ![net1_1](../assets/figures/snaqplot_net1_1.svg)
 
@@ -123,6 +125,7 @@ R"mtext"("hmax=2")           # add text annotation: title here
 plot(net3, :R, showGamma=true);
 R"mtext"("hmax=3")
 R"dev.off()"; # hide
+nothing # hide
 ```
 ![net23](../assets/figures/snaqplot_net23.svg)
 
@@ -272,6 +275,7 @@ R"svg(name('snaqplot_scores_heuristic.svg'), width=4, height=3)" # hide
 R"par"(mar=[2.5,2.5,.5,.5], mgp=[1.4,.4,0], tck=-0.02);  # hide
 R"plot"(scores, type="b", ylab="network score", xlab="hmax", col="blue");
 R"dev.off()"; # hide
+nothing # hide
 ```
 ![scores_heuristic](../assets/figures/snaqplot_scores_heuristic.svg)
 
@@ -306,6 +310,7 @@ R"svg(name('snaqplot_net1_2.svg'), width=4, height=3)" # starts image file
 R"par"(mar=[0,0,0,0]) # to reduce margins (no margins at all here)
 plot(net1, :R, showGamma=true, showEdgeNumber=true); # network is plotted & sent to file
 R"dev.off()"; # wrap up and save image file
+nothing # hide
 ```
 ![net1_2](../assets/figures/snaqplot_net1_2.svg)
 
@@ -322,6 +327,7 @@ R"svg(name('snaqplot_net1_3.svg'), width=4, height=3)" # hide
 R"par"(mar=[0,0,0,0]) # hide
 plot(net1, :R, showEdgeLength=true, minorHybridEdgeColor="tan")
 R"dev.off()"; # hide
+nothing # hide
 ```
 ![net1_3](../assets/figures/snaqplot_net1_3.svg)
 
@@ -342,6 +348,7 @@ R"par"(mar=[0,0,0,0]) # hide
 plot(net1,:R, tipOffset=0.5, showNodeNumber=true, edgeColor="tomato4",
      minorHybridEdgeColor="skyblue", majorHybridEdgeColor="tan");
 R"dev.off()"; # hide
+nothing # hide
 ```
 ![net1_4](../assets/figures/snaqplot_net1_4.svg)
 

src/auxiliary.jl

@@ -1210,44 +1210,59 @@ function switchHybridNode!(net::HybridNetwork, hybrid::Node, newHybrid::Node)
 end
 
 """
-`assignhybridnames!(net)`
+    assignhybridnames!(net)
 
-Assign names to hybrid nodes in the network `net`. Hybrid nodes with an empty `name` field ("")
-are modified with a name that does not conflict with other hybrid names in the network.
-The preferred name is "#H3" if the node number is 3 or -3, but an index other than 3 would be
-used if "#H3" were the name of another hybrid node already.
+Assign names to hybrid nodes in the network `net`.
+Hybrid nodes with an empty `name` field ("") are modified with a name that
+does not conflict with other hybrid names in the network. The preferred name
+is "H3" if the node number is 3 or -3, but an index other than 3 would be used
+if "H3" were the name of another node already.
 
 If two hybrid nodes have non-empty and equal names, the name of one of them is changed and
 re-assigned as described above (with a warning).
 """
 function assignhybridnames!(net::HybridNetwork)
-    hybnum = Int[]  # indices 'i' in hybrid names: #Hi
+    rx = r"^H(\d+)$"
+    # prep: collect indices 'i' of any tree nodes named like Hi
+    trenum = Int[]  # indices 'i' in tree node name, in case some are named Hi
+    for n in net.node
+        !n.hybrid || continue # do nothing if node n is hybrid
+        m = match(rx, n.name)
+        m == nothing || push!(trenum, parse(Int, m[1]))
+    end
     # first: go through *all* existing non-empty names
+    hybnum = Int[]  # indices 'i' in hybrid names: Hi
     for ih in 1:length(net.hybrid)
-        lab = net.hybrid[ih].name
+        hnode = net.hybrid[ih]
+        lab = hnode.name
         lab != "" || continue # do nothing if label is missing
         jh = findfirst(isequal(lab), [net.hybrid[j].name for j in 1:ih-1])
         if jh !== nothing # set repeated names to ""
-            @warn "hybrid nodes $(net.hybrid[ih].number) and $(net.hybrid[jh].number) have the same label: $lab. Will change the name of the former."
-            net.hybrid[ih].name = ""
-        else
-            m = match(r"^#H(\d+)$", lab)
-            if m != nothing # make full list of existing indices "i" in #Hi
-                push!(hybnum, parse(Int, m[1]))
+            @warn "hybrid nodes $(hnode.number) and $(net.hybrid[jh].number) have the same label: $lab. Will change the name of the former."
+            hnode.name = ""
+        else # fill in list of existing indices "i" in Hi
+            m = match(rx, lab)
+            m != nothing || continue # skip the rest if name is not of the form Hi
+            ind = parse(Int, m[1])
+            if ind in trenum
+                @warn "hybrid node $(hnode.number) had same label as a tree node: H$ind. Will change hybrid name."
+                hnode.name = ""
+            else
+                push!(hybnum, ind)
             end
         end
     end
-    # second: assign empty names to #Hi
+    # second: assign empty names to "Hi" for some i
     hnext = 1
     for ih in 1:length(net.hybrid)
         net.hybrid[ih].name == "" || continue # do nothing if non-empty label
         hnum = abs(net.hybrid[ih].number)
-        while in(hnum, hybnum)
+        while hnum in hybnum || hnum in trenum
             hnum = hnext  # not efficient, but rare
             hnext += 1    # and okay on small networks
         end
         push!(hybnum, hnum)
-        net.hybrid[ih].name = "#H$hnum"
+        net.hybrid[ih].name = "H$hnum"
     end
 end
 

src/readwrite.jl

@@ -47,33 +47,34 @@ end
 # allows names with letters and numbers: treats numbers as strings
 # it also reads # as part of the name and returns pound=true
 # it returns the node name as string as well to check if it exists already (as hybrid)
-function readNum(s::IO, c::Char, net::HybridNetwork, numLeft::Array{Int,1})
+function readnodename(s::IO, c::Char, net::HybridNetwork, numLeft::Array{Int,1})
     if !isValidSymbol(c)
         a = read(s, String);
-        error("Expected digit, alphanum or # at the start of taxon name, but received $(c). remaining is $(a).");
+        error("Expected digit, alphanum or # at the start of taxon name, but received $(c). remaining: $(a).");
     end
-    pound = (c == '#') ? 1 : 0
-    name = string(readskip!(s)) # reads in c from s, convert to string
-    c = peekskip(s)
+    pound = 0
+    name = ""
     while isValidSymbol(c)
-        d = readskip!(s)
-        name = string(name,d)
-        if d == '#'
-            pound += 1;
-            c = peekskip(s);
+        readskip!(s) # advance s past c (discard c, already read)
+        if c == '#'
+            pound += 1
+            c = readskip!(s) # read the character after the #
             if !isletter(c)
                 a = read(s, String);
-                error("Expected name after # but received $(c) in left parenthesis $(numLeft[1]-1). Remaining is $(a).")
+                error("Expected name after # but received $(c) in left parenthesis $(numLeft[1]-1). remaining: $(a).")
             end
             if c != 'L' && c != 'H' && c != 'R'
                 @warn "Expected H, R or LGT after # but received $(c) in left parenthesis $(numLeft[1]-1)."
             end
+            name = c * name # put the H, R or LGT first
+        else
+            name *= c # original c: put it last
         end
         c = peekskip(s);
     end
     if pound >1
         a = read(s, String);
-        error("strange node name with $(pound) # signs. remaining is $(a).")
+        error("strange node name with $(pound) # signs: $name. remaining: $(a).")
     end
     return size(net.names,1)+1, name, pound == 1
 end
@@ -397,13 +398,13 @@ function readSubtree!(s::IO, parent::Node, numLeft::Array{Int,1}, net::HybridNet
         c = peekskip(s);
         if isValidSymbol(c) # internal node has name
             hasname = true;
-            num, name, pound = readNum(s, c, net, numLeft);
+            num, name, pound = readnodename(s, c, net, numLeft);
             n.number = num; # n was given <0 number by parseRemainingSubtree!, now >0
             c = peekskip(s);
         end
     else # leaf, it should have a name
         hasname = true;
-        num, name, pound = readNum(s, c, net, numLeft)
+        num, name, pound = readnodename(s, c, net, numLeft)
         n = Node(num, true); # positive node number to leaves in the newick-tree description
         @debug "creating node $(n.number)"
     end
@@ -492,13 +493,13 @@ function readTopology(s::IO,verbose::Bool)
             elseif c == ')'
                 c = peekskip(s);
                 if isValidSymbol(c) # the root has a name
-                    num, name, pound = readNum(s, c, net, numLeft);
+                    num, name, pound = readnodename(s, c, net, numLeft);
                     n.name = name
                     # log warning or error if pound > 0?
                     c = peekskip(s);
                 end
                 if(c == ':') # skip information on the root edge, if it exists
-                    @warn "root edge ignored"
+                    # @warn "root edge ignored"
                     while c != ';'
                         c = readskip!(s)
                     end
@@ -932,18 +933,20 @@ Use `internallabel=false` to suppress the labels of internal nodes.
 """
 function writeSubTree!(s::IO, net::HybridNetwork, di::Bool, namelabel::Bool,
                        roundBL::Bool, digits::Integer, internallabel::Bool)
-    if net.numNodes == 1
-        print(s, (namelabel ? net.node[net.root].name : string(net.node[net.root].number)))
-    elseif net.numNodes > 1
+    rootnode = net.node[net.root]
+    if net.numNodes > 1
         print(s,"(")
-        degree = length(net.node[net.root].edge)
-        for e in net.node[net.root].edge
-            writeSubTree!(s,getOtherNode(e,net.node[net.root]),e,di,namelabel,roundBL,digits,internallabel)
+        degree = length(rootnode.edge)
+        for e in rootnode.edge
+            writeSubTree!(s,getOtherNode(e,rootnode),e,di,namelabel,roundBL,digits,internallabel)
             degree -= 1
             degree == 0 || print(s,",")
         end
         print(s,")")
     end
+    if internallabel || net.numNodes == 1
+        print(s, (namelabel ? rootnode.name : rootnode.number))
+    end
     print(s,";")
     return nothing
 end
@@ -956,7 +959,7 @@ end
 Write the extended newick format of the sub-network rooted at
 `node` and assuming that `edge` is a parent of `node`.
 
-If `parent` is `nothing`, the edge attribute `isChild1` is used
+If the parent `edge` is `nothing`, the edge attribute `isChild1` is used
 and assumed to be correct to write the subtree rooted at `node`.
 This is useful to write a subtree starting at a non-root node.
 Example:
@@ -997,9 +1000,10 @@ function writeSubTree!(s::IO, n::Node, parent::Union{Edge,Nothing},
     end
     # node label:
     if parent != nothing && parent.hybrid
-        print(s, (namelabel ? n.name : string("#H",n.number)))
+        print(s, "#")
+        print(s, (namelabel ? n.name : string("H", n.number)))
         n.name != "" || parent.isMajor || @warn "hybrid node $(n.number) has no name"
-    elseif (n.leaf)
+    elseif internallabel || n.leaf
         print(s, (namelabel ? n.name : n.number))
     end
     # branch lengths and γ, if available:
@@ -1388,7 +1392,6 @@ function writeTopology(net::HybridNetwork, s::IO,
     # finally, write parenthetical format
     writeSubTree!(s,net,di,true,round,digits,internallabel)
     # namelabel = true: to print leaf & node names (labels), not numbers
-    ## printID = false: print all branch lengths, not just identifiable ones
 end
 
 

test/test_compareNetworks.jl

@@ -14,16 +14,16 @@ global net, tree
 netstr = "(((A:4.0,(B:1.0)#H1:1.1::0.9):0.5,(C:0.6,#H1:1.0::0.1):1.0):3.0,D:5.0);"
 net = readTopology(netstr)
 @test_logs PhyloNetworks.deleteHybridEdge!(net, net.edge[6], true);
-@test writeTopology(net) == "(((A:4.0,(B:1.0):1.1):0.5,(C:0.6):1.0):3.0,D:5.0);"
+@test writeTopology(net) == "(((A:4.0,(B:1.0)H1:1.1):0.5,(C:0.6):1.0):3.0,D:5.0);"
 @test net.edge[3].gamma == 0.9
-@test net.node[3].name == "#H1"
+@test net.node[3].name == "H1"
 net = readTopology(netstr)
 @test_logs PhyloNetworks.deleteHybridEdge!(net, net.edge[3], true);
-@test writeTopology(net) == "(((A:4.0):0.5,(C:0.6,(B:1.0):1.0):1.0):3.0,D:5.0);"
+@test writeTopology(net) == "(((A:4.0):0.5,(C:0.6,(B:1.0)H1:1.0):1.0):3.0,D:5.0);"
 @test  net.edge[5].gamma == 0.1
 @test !net.edge[5].hybrid
 @test  net.edge[5].isMajor
-@test  net.node[3].name == "#H1"
+@test  net.node[3].name == "H1"
 @test !net.node[3].hybrid
 
 # example of network with one hybrid edge connected to the root
@@ -42,7 +42,7 @@ net = readTopology(netstr);
 @test writeTopology(net) == "(Adif:1.0,(Aech:0.122,(Asub:1.0,Agem:1.0):10.0):10.0);"
 net = readTopology(netstr);
 @test_logs PhyloNetworks.deleteHybridEdge!(net, net.edge[9], true);
-@test writeTopology(net) == "((Adif:1.0,(Aech:0.122,((Asub:1.0,Agem:1.0):0.0):10.0):10.0):1.614);"
+@test writeTopology(net) == "((Adif:1.0,(Aech:0.122,((Asub:1.0,Agem:1.0):0.0)H6:10.0):10.0):1.614);"
 
 if doalltests
 net=readTopology("(4,((1,(2)#H7:::0.864):2.069,(6,5):3.423):0.265,(3,#H7:::0.1361111):10.0);");
@@ -213,8 +213,8 @@ a = displayedTrees(net5, 0.1);
 
 net = readTopology("(((A:4.0,(B:1.0)#H1:1.1::0.9):0.5,(C:0.6,#H1:1.0::0.1):1.0):3.0,D:5.0);")
 trees = (@test_logs displayedTrees(net,0.0; keepNodes=true));
-@test writeTopology(trees[1])=="(((A:4.0,(B:1.0):1.1):0.5,(C:0.6):1.0):3.0,D:5.0);"
-@test writeTopology(trees[2])=="(((A:4.0):0.5,(C:0.6,(B:1.0):1.0):1.0):3.0,D:5.0);"
+@test writeTopology(trees[1])=="(((A:4.0,(B:1.0)H1:1.1):0.5,(C:0.6):1.0):3.0,D:5.0);"
+@test writeTopology(trees[2])=="(((A:4.0):0.5,(C:0.6,(B:1.0)H1:1.0):1.0):3.0,D:5.0);"
 @test PhyloNetworks.inheritanceWeight.(trees) ≈ [log(0.9), log(0.1)]
 
 end # of testset, displayedNetworks! & displayedTrees
@@ -226,9 +226,9 @@ net5 = readTopology("(A:1.0,((B:1.1,#H1:0.2::0.2):1.2,(((C:0.52,(E:0.5)#H2:0.02:
 @test_logs displayedNetworkAt!(net5, net5.hybrid[1]);
 @test writeTopology(net5) == "(A:1.0,((((C:0.52,(E:0.5)#H2:0.02::0.7):0.6,(#H2:0.01::0.3,F:0.7):0.8):0.9,D:1.1):1.3,B:2.3):0.7);"
 net = readTopology("((((B)#H1)#H2,((D,C,#H2)S1,(#H1,A)S2)S3)S4);") # missing γ's, level 2
-@test writeTopology(majorTree(net)) == "(((D,C),A),B);"
+@test writeTopology(majorTree(net)) == "(((D,C)S1,A)S3,B)S4;"
 setGamma!(net.edge[8], 0.8)
-@test writeTopology(majorTree(net)) == "((D,C),(A,B));"
+@test writeTopology(majorTree(net)) == "((D,C)S1,(A,B)S2)S3;"
 
 end # of testset, majorTree & displayedNetworkAt!
 

test/test_manipulateNet.jl

@@ -28,8 +28,8 @@ close(s); s = IOBuffer()
 @test String(take!(s)) == """
 node leaf  hybrid hasHybEdge name       inCycle edges'numbers
 1    true  false  false      B          -1      1   
-2    false true   true       #H1        -1      1    2    8   
-3    false true   true       #H2        -1      2    3    6   
+2    false true   true       H1         -1      1    2    8   
+3    false true   true       H2         -1      2    3    6   
 4    true  false  false      D          -1      4   
 5    true  false  false      C          -1      5   
 6    false false  true       S1         -1      4    5    6    7   
@@ -131,7 +131,7 @@ net.root=15; # node number -5 (clau: previously -4)
 @test_logs rootatnode!(net, -13); # or: rootatnode complained...
 @test_throws PhyloNetworks.RootMismatch rootatnode!(net, -5); # verbose = true this time
 # 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"; verbose=false); #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");

test/test_relaxed_reading.jl

@@ -22,7 +22,7 @@ global net
 	@test writeTopology(resultant) == writeTopology(expected)
 end
 @testset "isMajor and gamma consistency" begin
-	net = readTopology("((((B)#H1)#H2,((D,C,#H2:::0.8)S1,(#H1,A)S2)S3)S4);");
+	net = readTopology("((((B)#H1)#H2,((D,C,#H2:::0.8),(#H1,A))));");
 	@test writeTopology(net, round=true, digits=8) == "(#H2:::0.2,((D,C,((B)#H1)#H2:::0.8),(#H1,A)));"
 	net = readTopologyLevel1("(E,((B)#H1:::.5,((D,C),(#H1:::.5,A))));");
 	@test writeTopology(net) == "(D:1.0,C:1.0,((#H1:1.0::0.5,A:1.0):1.0,((B:1.0)#H1:1.0::0.5,E:1.0):1.0):1.0);"
@@ -32,12 +32,12 @@ end
 	redirect_stdout(originalstdout)
 end
 @testset "internal nodes" begin
-	net = (@test_logs (:warn, r"^root edge") readTopology("(a,b):0.5;"));
-	@test writeTopology(net) == "(a,b);"
-	@test writeTopology(readTopology("((a,(b)#H1)i1,(#H1,c)i2)r;")) == "((a,(b)#H1),(#H1,c));"
-	#@test writeTopology(readTopology("((a,(b)#H1)i1,(#H1,c))r;"), internallabel=true) ==
-	#	"((a,(b)#H1)i1,(#H1,c))r;"
-	@test_logs (:warn, r"^root edge") readTopology("((a,(b)#H1)i1,(#H1,c)i2)root:0.5;");
-	net = (@test_logs readTopology("(((a,(b)#H1)i1,(#H1,c)i2)root:0.5);"));
+	@test writeTopology(readTopology("(a,b):0.5;")) == "(a,b);"
+	@test writeTopology(readTopology("((a,(b)#H1)i1,(#H1,c))r;")) == "((a,(b)#H1)i1,(#H1,c))r;"
+	@test writeTopology(readTopology("((a,(b)#H1)i1,(#H1,c))r;"), internallabel=false) == "((a,(b)#H1),(#H1,c));"
+	@test_logs readTopology("((a,(b)#H1)i1,(#H1,c)i2)root:0.5;");
+	@test_logs readTopology("(((a,(b)#H1)i1,(#H1,c)i2)root:0.5);"); # root edge was deleted
+	# writeTopology(net) == "((a,(b)#H1)i1,(#H1,c)i2);"
+	# readTopology("((((a,(b)#H1)i1,(#H1,c)i2)root:0.5));"); still has 1 (of the 2) root edges
 end
 end

test/test_traitLikDiscrete.jl

@@ -71,7 +71,7 @@ a,b = randomTrait(m1, net2; keepInternal=false)
 Random.seed!(12345);
 a,b = randomTrait(m1, net2; keepInternal=true)
 @test a == [1  2  1  1  1  1  1  1  1]
-@test b == ["-2", "D", "-3", "-6", "C", "-4", "#H1", "B", "A"]
+@test b == ["-2", "D", "-3", "-6", "C", "-4", "H1", "B", "A"]
 if runall
     for e in net2.edge
         if e.hybrid 
@@ -262,7 +262,7 @@ fit1.model.rate[2] = 0.34981109618902395;
 asr = ancestralStateReconstruction(fit1)
 @test names(asr) == [:nodenumber, :nodelabel, :lo, :hi]
 @test asr[:nodenumber] == collect(1:9)
-@test asr[:nodelabel] == ["A","B","C","D","5","6","7","8","#H1"]
+@test asr[:nodelabel] == ["A","B","C","D","5","6","7","8","H1"]
 @test asr[:lo] ≈ [1.,1.,0.,0., 0.28602239466671175, 0.31945742289603263,
     0.16855042517785512, 0.7673588716207436, 0.7827758475866091] atol=1e-5
 @test asr[:hi] ≈ [0.,0.,1.,1.,0.713977605333288, 0.6805425771039674,