Multiple tree viz

.gitignore

@@ -3,4 +3,5 @@
 *.jl.mem
 /Manifest.toml
 /docs/build/
+/docs/src/generated/
 .DS_Store

docs/Project.toml

@@ -1,7 +1,10 @@
 [deps]
+Cairo = "159f3aea-2a34-519c-b102-8c37f9878175"
 Compose = "a81c6b42-2e10-5240-aca2-a61377ecd94b"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 FASTX = "c2308a5c-f048-11e8-3e8a-31650f418d12"
+Fontconfig = "186bb1d3-e1f7-5a2c-a377-96d770f13627"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MolecularEvolution = "9f975960-e239-4209-8aa0-3d3ad5a82892"
 Phylo = "aea672f4-3940-5932-aa44-993d1c3ff149"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"

docs/make.jl

@@ -1,10 +1,25 @@
 using MolecularEvolution
-using Documenter
+using Documenter, Literate
 using Phylo
 using Plots
-using Compose
+using Compose, Cairo, Fontconfig
 using FASTX
 
+LITERATE_INPUT = joinpath(@__DIR__, "..", "examples")
+LITERATE_OUTPUT = OUTPUT = joinpath(@__DIR__, "src/generated")
+
+for (root, _, files) ∈ walkdir(LITERATE_INPUT), file ∈ files
+    @show root, files
+    # ignore non julia files
+    splitext(file)[2] == ".jl" || continue
+    # full path to a literate script
+    ipath = joinpath(root, file)
+    # generated output path
+    opath = splitdir(replace(ipath, LITERATE_INPUT=>LITERATE_OUTPUT))[1]
+    # generate the markdown file calling Literate
+    Literate.markdown(ipath, opath)
+end
+
 DocMeta.setdocmeta!(
     MolecularEvolution,
     :DocTestSetup,
@@ -32,7 +47,7 @@ makedocs(;
         "simulation.md",
         "optimization.md",
         "ancestors.md",
-        "viz.md",
+        "generated/viz.md",
     ],
 )
 

examples/viz.jl

@@ -0,0 +1,103 @@
+# # Visualization
+
+# We offer two routes to visualization. The first is using our own plotting routines, built atop Compose.jl. The second converts our trees to Phylo.jl trees, and plots with their Plots.jl recipes. The Compose, Plots, and Phylo dependencies are optional.
+
+# ## Example 1
+
+using MolecularEvolution, Plots, Phylo
+
+#First simulate a tree, and then Brownian motion:
+tree = sim_tree(n = 20)
+internal_message_init!(tree, GaussianPartition())
+bm_model = BrownianMotion(0.0, 0.1)
+sample_down!(tree, bm_model)
+
+#We'll add the Gaussian means to the node_data dictionaries
+for n in getnodelist(tree)
+    n.node_data = Dict(["mu" => n.message[1].mean])
+end
+
+#Transducing the mol ev tree to a Phylo.jl tree
+phylo_tree = get_phylo_tree(tree)
+
+pl = plot(
+    phylo_tree,
+    showtips = true,
+    tipfont = 6,
+    marker_z = "mu",
+    markeralpha = 0.5,
+    line_z = "mu",
+    linecolor = :darkrainbow,
+    markersize = 4.0,
+    markerstrokewidth = 0,
+    margins = 1Plots.cm,
+    linewidth = 1.5,
+    markercolor = :darkrainbow,
+    size = (500, 500),
+)
+
+# We also offer `savefig_tweakSVG("simple_plot_example.svg", pl)` for some post-processing tricks that improve the exported trees, like rounding line caps, and `values_from_phylo_tree(phylo_tree,"mu")` which can extract stored quantities in the right order for passing into eg. `markersize` options when plotting.
+
+# For a more comprehensive list of things you can do with Phylo.jl plots, please see [their documentation](https://docs.ecojulia.org/Phylo.jl/stable/man/plotting/).
+
+# ## Drawing trees with `Compose.jl`.
+
+# The `Compose.jl` in-house tree drawing offers extensive flexibility. Here is an example that plots a pie chart representing the marginal probability of each of the 4 possible nucleotides on all nodes on the tree:
+
+using MolecularEvolution, Compose
+
+tree = sim_tree(40, 1000.0, 0.005, mutation_rate = 0.001)
+model = DiagonalizedCTMC(reversibleQ(ones(6), ones(4) ./ 4))
+internal_message_init!(tree, NucleotidePartition(ones(4) ./ 4, 1))
+sample_down!(tree, model)
+d = marginal_state_dict(tree, model);
+#-
+compose_dict = Dict()
+for n in getnodelist(tree)
+    compose_dict[n] =
+        (x, y) -> pie_chart(x, y, d[n][1].state[:, 1], size = 0.02, opacity = 0.75)
+end
+img = tree_draw(tree,draw_labels = false, line_width = 0.5mm, compose_dict = compose_dict)
+
+# This can then be exported with:
+
+savefig_tweakSVG("piechart_tree.svg",img);
+
+# ## Multiple trees
+# Doesn't require `Phylo.jl`. Query trees can be plotted against a reference tree with [`plot_multiple_trees`](@ref). This can be useful, for instance, when we've sampled trees with [`metropolis_sample`](@ref).
+using MolecularEvolution, Plots
+
+tree = sim_tree(10, 1, 1)
+nodelist = getnodelist(tree); mean = sum([n.branchlength for n in nodelist]) / length(nodelist)
+rparams(n::Int) = MolecularEvolution.sum2one(rand(n))
+model = DiagonalizedCTMC(reversibleQ(ones(6) ./ (6 * mean), rparams(4)))
+internal_message_init!(tree, NucleotidePartition(ones(4) ./ 4, 100))
+sample_down!(tree, model)
+@time trees, LLs = metropolis_sample(tree, [model], 300, collect_LLs=true);
+reference = trees[argmax(LLs)];
+# We'll use the maximum a posteriori tree as reference
+plot_multiple_trees(trees, reference)
+# We can pass in a weight function to fit query trees against `reference` in a weighted least squares fashion with a location and scale parameter.
+#=
+!!! note
+    If we don't want to scale the query trees, we must disable it with `opt_scale = false`.
+=#
+plot_multiple_trees(
+    trees,
+    reference,
+    y_jitter = 0.05,
+    weight_fn = n::FelNode ->
+       ifelse(MolecularEvolution.isroot(n) || isleafnode(n), 1.0, 0.0)
+)
+#jl # We only prefix with `MolecularEvolution` for `isroot`, since we have `Phylo.jl` in our namespace
+#=
+## Functions
+
+```@docs
+get_phylo_tree
+values_from_phylo_tree
+savefig_tweakSVG
+tree_draw
+plot_multiple_trees
+```
+=#

src/MolecularEvolution.jl

@@ -59,10 +59,11 @@ function __init__()
     end
 
     #Have Phylo and Plots installed if you want to draw trees with Phylo
-    @require Phylo = "aea672f4-3940-5932-aa44-993d1c3ff149" begin
-        @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
+    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
+        import .Plots
+        include("viz/multiple_trees.jl") #We're using "Plots." quite often so we might want to "using .Plots"
+        @require Phylo = "aea672f4-3940-5932-aa44-993d1c3ff149" begin
             import .Phylo
-            import .Plots
             include("viz/phylo_glue.jl")
         end
     end

src/utils/tree_hash.jl

@@ -7,13 +7,16 @@
 
 tuple_sort(t::Tuple{UInt64,UInt64}) = ifelse(t[1] < t[2], t, (t[2], t[1]))
 
-function node_hash_split(node, hash_container, node_container, name2hash)
+function node_hash_split(node, hash_container, node_container, name2hash; push_leaves = false)
     if isleafnode(node)
-        #consider if we want the children pairs too, and push in here too
+        if push_leaves
+            push!(hash_container, name2hash[node.name])
+            push!(node_container, node)
+        end
         return name2hash[node.name]
     else
         child_hashes = [
-            node_hash_split(nc, hash_container, node_container, name2hash) for
+            node_hash_split(nc, hash_container, node_container, name2hash, push_leaves = push_leaves) for
             nc in node.children
         ]
         #merge child hashes with xor as we go up the tree
@@ -27,15 +30,15 @@ function node_hash_split(node, hash_container, node_container, name2hash)
     end
 end
 
-function get_node_hashes(newt)
+function get_node_hashes(newt; push_leaves = false)
     leafnames = [n.name for n in getleaflist(newt)]
     leafhashes = hash.(leafnames)
     all_names_hash = xor(leafhashes...)
     name2hash = Dict(zip(leafnames, leafhashes))
     hash_container = UInt64[]
     node_container = FelNode[]
     #This puts things in the containers
-    node_hash_split(newt, hash_container, node_container, name2hash)
+    node_hash_split(newt, hash_container, node_container, name2hash, push_leaves = push_leaves)
     #This makes a hash that matches everything except the given node
     other_hash = xor.(hash_container, all_names_hash)
     #Sort these, to make comparisons order invariant, which makes the comparison rooting invariant
@@ -53,3 +56,31 @@ function tree_diff(query, reference)
     changed_nodes = newt_nc[[!(n in hashset) for n in newt_hc]]
     return changed_nodes
 end
+
+export tree_match_pairs
+function tree_match_pairs(query, reference; push_leaves = false)
+    newt_hc, newt_nc = get_node_hashes(query, push_leaves = push_leaves)
+    n_hc, n_nc = get_node_hashes(reference, push_leaves = push_leaves)
+    newt_hash2node = Dict(zip(newt_hc, newt_nc))
+    n_hash2node = Dict(zip(n_hc, n_nc))
+    return map(h -> (newt_hash2node[h], n_hash2node[h]), filter(x -> haskey(n_hash2node, x), newt_hc))
+end
+#=
+function tree_comp(query, reference, condition)
+    newt_hc, newt_nc = get_node_hashes(query)
+    n_hc, n_nc = get_node_hashes(reference)
+    changed_nodes = filter(condition(n_hc), newt_hc)
+    return changed_nodes
+end
+
+export tree_diff
+tree_diff(query, reference) = tree_comp(query, reference, !in ∘ Set)
+
+export tree_match
+function tree_match(query, reference)
+    filter(newt_hc, n_hc) = map(in(Set(n_hc)), newt_hc)
+    return tree_comp(query, reference, filter)
+end
+
+tree_match(query, reference) = tree_comp(query, reference, in ∘ Set)
+=#