utilities to bind phylogenetic networks, aka add or graft a subnetwork onto a larger network

addsubnetwork.jl

#=
to use the functions in this file, read it within julia with:
include("path/to/addsubnetwork.jl");
but replace 'path/to/' by the actual path to this file on your machine.
The main functions defined in this file are:

addnetwork_at!
addnetwork_atleaf!

=#

try # install package PhyloNetworks, if not already installed, and load it
    using PhyloNetworks
catch
    @info "need PhyloNetworks package: will add & use it"
    import Pkg; Pkg.add("PhyloNetworks"); using PhyloNetworks
end

"""
    addnetwork_at!

Add (graft or bind) to network `net` the subnetwork `subnet`, attaching it
to `node`, which should be a node in `net`.
`net` is modified. `subnet` is not.
- If `edgelength` is not provided, then `node` is replaced by the root of `subnet`.
- If `edgelength` is provided, an edge of length `edgelength` is created to form
  the stem of the subnetwork, from `node` to the root of `subnet`.
  If `edgelength=-1.0`, the stem edge length is interpreted as missing.
  If `verbose` is true, then a warning will appear in case `edgelength` is
  negative other than -1.

output: crown node of subnetwork in full network `net`.

# examples

```julia
julia> newick = "((S2c:3,(S1:1)#H1),(#H1,S4));"; net = readTopology(newick);

julia> subnet2 = readTopology("(S2a,S2b);"); # we want this clade to be sister to S2c

julia> leaf_S2c = net.node[findfirst(n -> n.name == "S2c", net.node)];

julia> edge_to_S2c = getparentedge(leaf_S2c)
PhyloNetworks.EdgeT{PhyloNetworks.Node}:
 number:1
 length:3.0
 attached to 2 node(s) (parent first): -3 1

julia> newnode, newedge = PhyloNetworks.breakedge!(edge_to_S2c, net);

julia> (newedge.length, edge_to_S2c.length)
(1.5, 1.5)

julia> crown = addnetwork_at!(net, newnode, subnet2);

julia> writeTopology(net) # polytomy of three S2 taxa
(((S1:1.0)#H1,(S2c:1.5,S2a,S2b):1.5),(#H1,S4));

julia> net = readTopology(newick); leaf_S2c = net.node[1]; # re-do, but without a polytomy

julia> newnode, _ = PhyloNetworks.breakedge!(getparentedge(leaf_S2c), net);

julia> crown = addnetwork_at!(net, newnode, subnet2, 0.5);

julia> writeTopology(net) # resolved: S2c | S2a,S2b
"(((S1:1.0)#H1,(S2c:1.5,(S2a,S2b)_subnetwork_crown:0.5):1.5),(#H1,S4));"

```
"""
function addnetwork_at!(
    net::HybridNetwork,
    node::PhyloNetworks.Node,
    subnet::HybridNetwork;
)
    subnet = deepcopy(subnet) # to leave input subnetwork unchanged
    crown = subnet.node[subnet.root]    # extract the crown node
    deleteat!(subnet.node, subnet.root) # all other nodes will transfer
    # re-number subnetwork nodes and edges, to make them unique in full network
    nextnum = maximum(n.number for n in net.node) + 1
    for n in subnet.node
        n.number = nextnum
        nextnum += 1
    end
    nextnum = maximum(e.number for e in net.edge) + 1
    for e in subnet.edge
        e.number = nextnum
        nextnum +=1
    end
    # if node was a leaf originally, then won't be any more
    if node.leaf && !isempty(crown.edge)
        deleteat!(net.leaf, findfirst(isequal(node), net.leaf))
        node.leaf = false
        net.numTaxa -= 1
    end
    for ee in crown.edge # detacth each edge adjacent to crown, attach to `node` instead
        PhyloNetworks.removeNode!(crown, ee) 
        PhyloNetworks.setNode!(ee, node) # node comes 2nd in ee now, so isChild1 most likely wrong
        PhyloNetworks.setEdge!(node, ee)
    end
    # append to net.{node,edge} all subnetworks's nodes and edges, except for crown
    for n in subnet.node PhyloNetworks.pushNode!(net, n); end
    for e in subnet.edge PhyloNetworks.pushEdge!(net, e); end
    directEdges!(net) # correct isChild1
    # cleanup, to help garbage collection
    empty!(crown.edge);   empty!(subnet.node);   empty!(subnet.edge);
    empty!(subnet.names); empty!(subnet.hybrid); empty!(subnet.leaf)
    return node
end
function addnetwork_at!(
    net::HybridNetwork,
    node::PhyloNetworks.Node,
    subnet::HybridNetwork,
    edgelength::Float64;
    verbose::Bool=true,
)
    verbose && edgelength < 0 && edgelength != -1 &&
        @warn "stem of negative length: $edgelength"
    crownname = node.name * "_subnetwork_crown"
    newleaf = PhyloNetworks.addleaf!(net, node, crownname, edgelength)
    return addnetwork_at!(net, newleaf, subnet)
end

"""
    addnetwork_atleaf!

Replace a leaf node `n` in `net` by the subnetwork `subnet`.
`n` is `leaf` if it's a node.
If `leaf` is a string then `n` is the (first) leaf node named `leaf`.

As `leaf` becomes internal, and often times the subnetwork may have a taxon of
the same name, this node is renamed by appending "_subnetwork_crown" to its
original name.

`edgelength` controls how the length of the external (or pendent) edge to the
leaf is modified. After the leaf replacement, this edge is the subnetwork's stem.
- if `edgelength` is a number, then the edge to the original leaf is given
  length `edgelength` (with -1 interpreted as missing).
  If `verbose` is true, then a warning will appear in case `edgelength` is
  negative other than -1.
- if `edgelength=nothing`, then the edge to the original leaf is suppressed.
  This may result in a polytomy in case `subnet` does not already have a stem
  edge (has a bifurcation root instead)

Uses `addnetwork_at!`.

# examples

```julia
julia> newick = "((S2,(S1:1)#H1),(#H1,S4));"

julia> subnet1 = readTopology("(S1a,S1b);"); # we want to replace S1 by this clade

julia> net = readTopology(newick);

julia> addnetwork_atleaf!(net, "S1", subnet1, 0.3); writeTopology(net)
"((S2,((S1a,S1b)S1_subnetwork_crown:0.3)#H1),(#H1,S4));"

julia> net = readTopology(newick); # re-do, but keep original length for the stem

julia> crown = addnetwork_atleaf!(net, "S1", subnet1, nothing; crownname="");

julia> writeTopology(net)
"((S2,((S1a,S1b):1.0)#H1),(#H1,S4));"

julia> stemedge = getparentedge(crown)

julia> PhyloNetworks.shrinkedge!(net, stemedge); # shrink the stem, if desired

julia> writeTopology(net) # here we get a polytomy below H1
"((S2,(S1a,S1b)#H1),(#H1,S4));"
```

In this other example, we have 2 networks: one large, one small. We do this below:
1. check that these two networks share exactly 2 leaves, one to be interpreted
   as an outgroup, the other to be swapped by the small network in large network.
2. delete the outgroup from the small network
3. add the small network to the large network, at their shared leaf.

```julia
julia> largenet = readTopology("(((S2a:3,(S1:1)#H1),(#H1,S4)),(O1,O2));");

julia> smallnet = readTopology("((S2a,S2b),O1);");

julia> sharedtiplabels = intersect(tipLabels(largenet), tipLabels(smallnet))
2-element Vector{String}:
 "S2a"
 "O1"

julia> deleteleaf!(smallnet, "O1") # delete outgroup from small network

julia> sharedtiplabels = intersect(tipLabels(largenet), tipLabels(smallnet))
1-element Vector{String}:
 "S2a"

julia> if length(sharedtiplabels) != 1 error("oops, I wanted a single shared leaf"); end

julia> leafname = sharedtiplabels[1]

julia> addnetwork_atleaf!(largenet, leafname, smallnet, nothing; crownname="smallnet_crown");

julia> writeTopology(largenet)
"((((S2a,S2b)smallnet_crown:3.0,(S1:1.0)#H1),(#H1,S4)),(O1,O2));"

```
"""
function addnetwork_atleaf!(n::HybridNetwork, leaf::String, args...; kwargs...)
    i = findfirst(x -> x.name == leaf, n.node)
    isnothing(i) && error("no node named $leaf to replace by a subnetwork")
    return addnetwork_atleaf!(n, n.node[i], args...; kwargs...)
end
function addnetwork_atleaf!(
    net::HybridNetwork,
    leaf::PhyloNetworks.Node,
    subnet::HybridNetwork,
    edgelength::Union{Number,Nothing};
    verbose::Bool=true,
    crownname = leaf.name * "_subnetwork_crown",
)
    leaf.hybrid &&
        error("node $(leaf.name) number $(leaf.number) is a hybrid. use addnetwork_at! instead?")
    leaf.leaf ||
        @error("node $(leaf.name) number $(leaf.number) is not a leaf")
    pe = getparentedge(leaf)
    leaf.name = crownname
    if !isnothing(edgelength)
        verbose && edgelength < 0 && edgelength != -1 &&
            @warn "stem of negative length: $edgelength"
        pe.length = edgelength
    end
    return addnetwork_at!(net, leaf, subnet)
end

println("""
    Done reading file. Function "addnetwork_atleaf!" requires PhyloNetworks v0.16.4.
    To get help, type  ?  then   addnetwork_atleaf!""")