text:
# functions to read/write networks topologies

# peek the next non-white-space char
# removes white spaces from the IOStream/IOBuffer
# see skipchars(predicate, io::IO; linecomment=nothing) in io.jl
# https://github.com/JuliaLang/julia/blob/3b02991983dd47313776091720871201f75f644a/base/io.jl#L971
ncodeunits(c::Char) = write(devnull, c) # for Julia v0.7. Remove for Julia v1.0.3+
function peekskip(io::IO, linecomment=nothing)
    c = missing
    while !eof(io)
        c = read(io, Char)
        if c === linecomment
            readline(io)
        elseif !isspace(c)
            skip(io, -ncodeunits(c))
            break
        end
    end
    return c # skipchar returns io instead
end

# aux function to read the next non-white-symbol char in s, advances s
# input: s IOStream/IOBuffer
function readskip!(io::IO)
    c = missing
    while !eof(io)
        c = read(io, Char)
        if !isspace(c)
            break
        end
    end
    return c
end

# advance stream in readSubtree
function advance!(s::IO, numLeft::Array{Int,1})
    c = readskip!(s)
    if ismissing(c)
        error("Tree ends prematurely while reading subtree after left parenthesis $(numLeft[1]-1).")
    end
    return c
end

# auxiliary function to read a taxon name.
# 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 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: $(a).");
    end
    pound = 0
    name = ""
    while isValidSymbol(c)
        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: $(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: $name. remaining: $(a).")
    end
    return size(net.names,1)+1, name, pound == 1
end

# aux function to read floats like length or gamma values, to be read after a colon
function readFloat(s::IO, c::Char)
    if !(isdigit(c) || c in ['.','e','-','E'])
        a = read(s, String);
        error("Expected float digit after ':' but found $(c). remaining is $(a).");
    end
    num = ""
    while isdigit(c) || c in ['.','e','-', 'E']
        d = read(s, Char) # reads c and advances IO
        num = string(num,d);
        c = peekskip(s);
    end
    f = 0.0
    try
        f = parse(Float64, num)
    catch
        error("problem with number read $(num), not a float number")
    end
    return f
end

# aux function to identify if a symbol in taxon name is valid
# allowed: letters, numbers, underscores, # (for hybrid name)
# NOT allowed: white space () [] : ; ' ,
# according to richnewick.pdf
# actually: ' is allowed. coded weirdly imo!
function isValidSymbol(c::Char)
    return isletter(c) || isnumeric(c) || c=='_' || c=='#' ||
      (!isspace(c) && c != '(' && c != ')' && c != '[' && c != ']' && c != ':' && c != ';' && c != ',')
end

"""
    parseRemainingSubtree!(s::IO, numLeft, net, hybrids)

Create internal node. Helper for [`readSubtree!`](@ref),
which creates the parent edge of the node created by `parseRemainingSubtree!`:
`readSubtree!` calls `parseRemainingSubtree!`, and vice versa.
Called once a `(` has been read in a tree topology and reads until the corresponding `)` has been found.
This function performs the recursive step for `readSubtree!`.
Advances `s` past the subtree, adds discovered nodes and edges to `net`, and `hybrids`.

Does *not* read the node name and the edge information of the subtree root:
this is done by [`readSubtree!`](@ref)
"""
@inline function parseRemainingSubtree!(s::IO, numLeft::Array{Int,1}, net::HybridNetwork, hybrids::Vector{String})
    numLeft[1] += 1
    DEBUGC && @debug "" numLeft
    n = Node(-1*numLeft[1],false);
    @debug "creating node $(n.number)"
    keepon = true;
    c = readskip!(s)
    while (keepon)
        bl = readSubtree!(s,n,numLeft,net,hybrids)
        c = advance!(s,numLeft)
        if c == ')'
            keepon = false
        elseif c != ','
            a = read(s, String);
            error("Expected right parenthesis after left parenthesis $(numLeft[1]-1) but read $(c). The remainder of line is $(a).")
        end
    end
    return n
end

"""
    parseHybridNode!(node, parentNode, hybridName, net, hybrids)

Helper function for `readSubtree!`. Create the parent edge for `node`.
Return this edge, and the hybrid node retained (`node` or its clone in the newick string).
Insert new edge and appropriate node into `net` and `hybrids` accordingly.
Handles any type of given hybrid node.
Called after a `#` has been found in a tree topology.
"""
@inline function parseHybridNode!(n::Node, parent::Node, name::String, net::HybridNetwork, hybrids::Vector{String})
    @debug "found pound in $(name)"
    n.hybrid = true;
    DEBUGC && @debug "got hybrid $(name)"
    DEBUGC && @debug "hybrids list has length $(length(hybrids))"
    ind = findfirst(isequal(name), hybrids) # index of 'name' in the list 'hybrid'. nothing if not found
    e = Edge(net.numEdges+1) # isMajor = true by default
    if n.leaf e.isMajor = false; end
    e.hybrid = true
    e.gamma = -1.0
    if ind !== nothing # the hybrid name was seen before
        @debug "$(name) was found in hybrids list"
        ni = findfirst(isequal(name), [no.name for no in net.node])
        ni !== nothing || error("hybrid name $name was supposed to be in the network, but not found")
        other = net.node[ni]
        @debug "other is $(other.number)"
        DEBUGC && @debug "other is leaf? $(other.leaf), n is leaf? $(n.leaf)"
        if !n.leaf && !other.leaf
            error("both hybrid nodes are internal nodes: successors of the hybrid node must only be included in the node list of a single occurrence of the hybrid node.")
        elseif n.leaf
            @debug "n is leaf"
            @debug "creating hybrid edge $(e.number) attached to other $(other.number) and parent $(parent.number)"
            pushEdge!(net,e);
            setNode!(e,[other,parent]); # isChild1 = true by default constructor
            setEdge!(other,e);
            setEdge!(parent,e);
            n = other # original 'n' dropped, 'other' retained: 'n' output to modify 'n' outside
            @debug "e $(e.number )istIdentifiable? $(e.istIdentifiable)"
        else # !n.leaf : delete 'other' from the network
            @debug "n is not leaf, other is leaf"
            size(other.edge,1) == 1 || # other should be a leaf
               error("strange: node $(other.number) is a leaf hybrid node. should have only 1 edge but has $(size(other.edge,1))")
            DEBUGC && @debug "other is $(other.number), n is $(n.number), edge of other is $(other.edge[1].number)"
            otheredge = other.edge[1];
            otherparent = getOtherNode(otheredge,other);
            @debug "otheredge is $(otheredge.number)"
            @debug "parent of other is $(otherparent.number)"
            removeNode!(other,otheredge);
            deleteNode!(net,other);
            setNode!(otheredge,n);
            setEdge!(n,otheredge);
            ## otheredge.istIdentifiable = true ## setNode should catch this, but when fixed, causes a lot of problems
            @debug "setting otheredge to n $(n.number)"
            @debug "creating hybrid edge $(e.number) between n $(n.number) and parent $(parent.number)"
            setNode!(e,[n,parent]);
            setEdge!(n,e);
            setEdge!(parent,e);
            pushNode!(net,n);
            pushEdge!(net,e);
            n.number = other.number; # modifies original negative node number, to positive node #
            n.name = other.name;
            @debug "edge $(e.number) istIdentifiable? $(e.istIdentifiable)"
            @debug "otheredge $(otheredge.number) istIdentifiable? $(otheredge.istIdentifiable)"
        end
    else # ind==nothing: hybrid name not seen before
        @debug "$(name) not found in hybrids list"
        @debug "$(name) is leaf? $(n.leaf)"
        n.hybrid = true;
        nam = string(name)
        push!(net.names, nam);
        n.name = nam;
        DEBUGC && @debug "put $(nam) in hybrids name list"
        push!(hybrids, nam);
        pushNode!(net,n);
        @debug "creating hybrid edge $(e.number)"
        pushEdge!(net,e);
        setNode!(e,[n,parent]);
        setEdge!(n,e);
        setEdge!(parent,e);
        @debug "edge $(e.number) istIdentifiable? $(e.istIdentifiable)"
    end
    e.containRoot = !e.hybrid # not good: but necessay for SNaQ functions
    return (e,n)
end

"""
    parseTreeNode!(node, parentNode, net)

Helper function for `readSubtree!`.
Insert the input tree node and associated edge (created here) into `net`.
"""
@inline function parseTreeNode!(n::Node, parent::Node, net::HybridNetwork)
    pushNode!(net,n);
    e = Edge(net.numEdges+1);
    pushEdge!(net,e);
    setNode!(e,[n,parent]);
    setEdge!(n,e);
    setEdge!(parent,e);
    return e
end

"""
    getdataValue!(s::IO, int, numLeft::Array{Int,1})

Helper function for `parseEdgeData!`.
Read a single floating point edge data value in a tree topology.
Return -1.0 if no value exists before the next colon, return the value as a float otherwise.
Modifies s by advancing past the next colon character.
Only call this function to read a value when you know a numerical value exists!
"""
@inline function getDataValue!(s::IO, call::Int, numLeft::Array{Int,1})
    errors = ["one colon read without double in left parenthesis $(numLeft[1]-1), ignored.",
              "second colon : read without any double in left parenthesis $(numLeft[1]-1), ignored.",
              "third colon : without gamma value after in $(numLeft[1]-1) left parenthesis, ignored"]
    c = peekskip(s)

    # Value is present
    if isdigit(c) || c == '.'
        val = readFloat(s, c)
        return val
    # No value
    elseif c == ':'
        return -1.0
    else
        @warn errors[call]
        return -1.0
    end
end

"""
    parseEdgeData!(s::IO, edge, node, numberOfLeftParentheses::Array{Int,1})

Helper function for readSubtree!, fixes a bug from using setGamma
Modifies `e` according to the specified edge length and gamma values in the tree topology.
Advances the stream `s` past any existing edge data.
Edges in a topology may optionally be followed by ":edgeLen:bootstrap:gamma"
where edgeLen, bootstrap, and gamma are decimal values.
"""
@inline function parseEdgeData!(s::IO, e::Edge, n::Node, numLeft::Array{Int,1})
    read(s, Char); # to read the first ":"
    e.length = getDataValue!(s, 1, numLeft)
    bootstrap = nothing;
    if peekskip(s) == ':'
        readskip!(s)
        bootstrap = getDataValue!(s, 2, numLeft)
    end
    # e.gamma = -1.0 by default when e is created by parseHybridNode!
    if peekskip(s) == ':'
        readskip!(s)
        e.gamma = getDataValue!(s, 3, numLeft)
    end
    if e.gamma != -1.0 && !e.hybrid && e.gamma != 1.0
        @warn "γ read for edge $(e.number) but it is not hybrid, so γ=$(gamma) ignored"
        e.gamma = 1.0
    end
end

"""
    synchronizePartnersData!(e::Edge, n::Node)

Synchronize γ and isMajor for edges `e` and its partner,
both hybrid edges with the same child `n`:

- if one γ is missing and the other is not: set the missing γ to 1 - the other
- γ's should sum up to 1.0
- update `isMajor` to match the γ information: the major edge is the one with γ > 0.5.

**Warnings**: does not check that `e` is a hybrid edge,
nor that `n` is the child of `e`.
"""
@inline function synchronizePartnersData!(e::Edge, n::Node)
    partners = Edge[] # The edges having n as a child, other than e
    for e2 in n.edge
        if e2.hybrid && e2!=e && n==getChild(e2)
            push!(partners, e2)
        end
    end
    numPartners = length(partners)
    if numPartners == 0
        return nothing
    end
    if numPartners > 1 # 3 or more hybrid parents
        error("γ value parsed but hybrid edge has $numPartners partners (should be 0 or 1)")
    end
    # numPartners == 1 then. partner edge has been read, may have gamma set
    partner = partners[1]
    if e.gamma == -1. && partner.gamma == -1. # exact -1.0 means missing
        return nothing
    end
    # γ non-missing for e and/or partner, then
    # update γ and isMajor of both edges, to be consistent with each other
    if e.gamma == -1.
        if partner.gamma < 0.0
            @warn "partners: $(e.number) with no γ, $(partner.number) with γ<0. will turn to 1 and 0"
            partner.gamma = 0.0
            e.gamma = 1.0
        elseif partner.gamma > 1.0
            @warn "partners: $(e.number) with no γ, $(partner.number) with γ>1. will turn to 0 and 1"
            partner.gamma = 1.0
            e.gamma = 0.0
        else # partner γ in [0,1]
            e.gamma = 1. - partner.gamma
        end
    elseif partner.gamma == -1.
        if e.gamma < 0.0
            @warn "partners: $(partner.number) with no γ, $(e.number) with γ<0. will turn to 1 and 0"
            e.gamma = 0.0
            partner.gamma = 1.0
        elseif e.gamma > 1.0
            @warn "partners: $(partner.number) with no γ, $(e.number) with γ>1. will turn to 0 and 1"
            e.gamma = 1.0
            partner.gamma = 0.0
        else # γ in [0,1]
            partner.gamma = 1. - e.gamma
        end
    else # both γ's are non-missing. won't check if negative
        gammasum = e.gamma + partner.gamma
        if !isapprox(gammasum, 1.0)
            e.gamma = e.gamma/gammasum
            partner.gamma = 1. - e.gamma
        end
    end
    # next: update isMajor, originally based on which node was a leaf and which was not
    # if both γ are 0.5: keep isMajor as is. Otherwise: γ's take precedence.
    if !isapprox(e.gamma, 0.5)
        emajor = e.gamma > 0.5
        if e.isMajor != emajor # leaf status was inconsistent with γ info
            e.isMajor = emajor
            partner.isMajor = !emajor
        end
    end
end

"""
    readSubtree!(s::IO, parentNode, numLeft, net, hybrids)

Recursive helper method for `readTopology`:
read a subtree from an extended Newick topology.
input `s`: IOStream/IOBuffer.

Reads additional info formatted as: `:length:bootstrap:gamma`.
Allows for name of internal nodes without # after closing parenthesis: (1,2)A.
Warning if hybrid edge without γ, or if γ (ignored) without hybrid edge
"""
function readSubtree!(s::IO, parent::Node, numLeft::Array{Int,1}, net::HybridNetwork, hybrids::Vector{String})
    c = peekskip(s)
    e = nothing;
    hasname = false; # to know if the current node has name
    pound = false;
    if c == '('
        # read the rest of the subtree (perform the recursive step!)
        n = parseRemainingSubtree!(s, numLeft, net, hybrids)
        c = peekskip(s);
        if isValidSymbol(c) # internal node has name
            hasname = true;
            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 = 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
    if pound # found pound sign in name
        # merge the 2 nodes corresponding the hybrid: make n the node that is retained
        e,n = parseHybridNode!(n, parent, name, net, hybrids) # makes hybrid node number >0
    else
        if hasname
            push!(net.names, name);
            n.name = name
        end
        e = parseTreeNode!(n, parent, net)
    end
    c = peekskip(s);
    e.length = -1.0
    if c == ':'
        parseEdgeData!(s, e, n, numLeft)
    end
    if e.hybrid
        # if hybrid edge: 'e' might have no info, but its partner may have had info
        synchronizePartnersData!(e, n) # update γ and isMajor of e and/or its partner
    end
    return true
end

# function to read topology from parenthetical format
# input: file name or tree in parenthetical format
# calls readTopology(s::IO)
# warning: crashes if file name starts with (
function readTopology(input::AbstractString,verbose::Bool)
    if(input[1] == '(') # input = parenthetical description
       s = IOBuffer(input)
    else # input = file name
        try
            s = open(input)
        catch
            error("Could not find or open $(input) file");
        end
       s = open(input)
    end
    net = readTopology(s,verbose)
    return net
end

"""
    readTopology(file name)
    readTopology(parenthetical description)

Read tree or network topology from parenthetical format (extended Newick).
If the root node has a single child: ignore (i.e. delete from the topology)
the root node and its child edge.

Input: text file or parenthetical format directly.
The file name may not start with a left parenthesis, otherwise the file
name itself would be interpreted as the parenthetical description.

A root edge, not enclosed within a pair a parentheses, is ignored.
If the root node has a single edge, this one edge is removed.
"""
readTopology(input::AbstractString) = readTopology(input,true)

function readTopology(s::IO,verbose::Bool)
    net = HybridNetwork()
    line = readuntil(s,";", keep=true);
    if(line[end] != ';')
        error("file does not end in ;")
    end
    seekstart(s)
    c = peekskip(s)
    numLeft = [1]; # made Array to make it mutable; start at 1 to avoid node -1 which breaks undirectedOtherNetworks
    hybrids = String[];
    if(c == '(')
        numLeft[1] += 1;
        #println(numLeft)
        n = Node(-1*numLeft[1],false);
        c = readskip!(s)
        b = false;
        while(c != ';')
            b |= readSubtree!(s,n,numLeft,net,hybrids)
            c = readskip!(s);
            if eof(s)
                error("Tree ended while reading in subtree beginning with left parenthesis number $(numLeft[1]-1).")
            elseif c == ','
                continue;
            elseif c == ')'
                c = peekskip(s);
                if isValidSymbol(c) # the root has a name
                    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"
                    while c != ';'
                        c = readskip!(s)
                    end
                end
            end
        end
        @debug "after readsubtree:"
        @debug begin printEdges(net); "printed edges" end
        # delete the root edge, if present
        if size(n.edge,1) == 1 # root node n has only one edge
            edge = n.edge[1]
            child = getOtherNode(edge,n);
            removeEdge!(child,edge);
            net.root = getIndex(child,net);
            deleteEdge!(net,edge);
        else
            pushNode!(net,n);
            net.root = getIndex(n,net);
        end
    else
		a = read(s, String)
        error("Expected beginning of tree with ( but received $(c) instead, rest is $(a)")
    end
    storeHybrids!(net)
    checkNumHybEdges!(net)
    directEdges!(net; checkMajor=true) # to update edges containRoot: true until hybrid, false below hybrid
    net.isRooted = true
    return net
end

readTopology(s::IO) = readTopology(s,true)

"""
    `checkNumHybEdges!(net)`

Check for consistency between hybrid-related attributes in the network:
- for each hybrid node: 2 or more hybrid edges
- exception: allows for a leaf to be attached to a single hybrid edge
- exactly 2 incoming parent hybrid edges
Run after `storeHybrids!`. See also `check2HybEdges`.
"""
function checkNumHybEdges!(net::HybridNetwork)
    if isTree(net) return nothing; end
    !isempty(net.hybrid) || error("net.hybrid should not be empty for this network")
    for n in net.hybrid
        hyb = sum([e.hybrid for e in n.edge]); # number of hybrid edges attached to node
        if hyb == 1
            if net.numHybrids == 1
                error("only one hybrid node $(n.number) named $(n.name) found with one hybrid edge attached")
            else
                error("hybrid node $(n.number) named $(n.name) has only one hybrid edge attached. there are $(net.numHybrids-1) other hybrids out there but this one remained unmatched")
            end
        elseif hyb == 0
            if length(n.edge) == 0
                error("strange hybrid node $(n.number) attached to 0 edges")
            elseif length(n.edge) == 1
                n.leaf || error("hybrid node $(n.number) has only one tree edge attached and it is not a leaf")
            elseif length(n.edge) >= 2
                @warn "hybrid node $(n.number) not connected to any hybrid edges. Now transformed to tree node"
                n.hybrid = false;
            end
        elseif hyb >=2 # check: exactly 2 incoming, no more.
            nhybparents = 0
            for e in n.edge
                if n == getChild(e)
                    if e.hybrid
                        nhybparents += 1
                    else @error "node $(n.number) has parent tree edge $(e.number): wrong isChild1 for this edge?"
                    end
                end
            end
            if nhybparents < 2
                error("hybrid node $(n.number) with $nhybparents = fewer than 2 hybrid parents")
            elseif nhybparents >2
                error("hybrid node $(n.number) with $nhybparents: we don't allow such polytomies")
            end
        end
    end
    return nothing
end

# aux function to send an error if the number of hybrid attached to every
# hybrid node is >2
function check2HybEdges(net::HybridNetwork)
    for n in net.hybrid
        hyb = sum([e.hybrid for e in n.edge]);
        if hyb > 2
            error("hybrid node $(n.number) has more than two hybrid edges attached to it: polytomy that cannot be resolved without intersecting cycles.")
        end
    end
end

# aux function to solve a polytomy
# warning: chooses one resolution at random
function solvePolytomyRecursive!(net::HybridNetwork, n::Node)
    if(size(n.edge,1) == 4)
        edge1 = n.edge[1];
        edge2 = n.edge[2];
        edge3 = n.edge[3];
        edge4 = n.edge[4];
        removeEdge!(n,edge3);
        removeEdge!(n,edge4);
        removeNode!(n,edge3);
        removeNode!(n,edge4);
        ednew = Edge(net.numEdges+1,0.0);
        max_node = maximum([e.number for e in net.node]);
        n1 = Node(max_node+1,false,false,[edge3,edge4,ednew]);
        setEdge!(n,ednew);
        setNode!(edge3,n1);
        setNode!(edge4,n1);
        setNode!(ednew,[n,n1]);
        pushNode!(net,n1);
        pushEdge!(net,ednew);
    else
        edge1 = n.edge[1];
        removeEdge!(n,edge1);
        solvePolytomyRecursive!(net,n);
        setEdge!(n,edge1);
    end
end

# function to solve a polytomy among tree edges recursively
function solvePolytomy!(net::HybridNetwork, n::Node)
    !n.hybrid || error("cannot solve polytomy in a hybrid node $(n.number).")
    while(size(n.edge,1) > 3)
        solvePolytomyRecursive!(net,n);
    end
end

# aux function to add a child to a leaf hybrid
function addChild!(net::HybridNetwork, n::Node)
    n.hybrid || error("cannot add child to tree node $(n.number).")
    ed1 = Edge(net.numEdges+1,0.0);
    n1 = Node(size(net.names,1)+1,true,false,[ed1]);
    setEdge!(n,ed1);
    setNode!(ed1,[n,n1]);
    pushNode!(net,n1);
    pushEdge!(net,ed1);
end
# aux function to expand the children of a hybrid node
function expandChild!(net::HybridNetwork, n::Node)
    if n.hybrid
        suma = count([!e.hybrid for e in n.edge]);
        #println("create edge $(net.numEdges+1)")
        ed1 = Edge(net.numEdges+1,0.0);
        n1 = Node(size(net.names,1)+1,false,false,[ed1]);
        #println("create node $(n1.number)")
        hyb = Edge[];
        for i in 1:size(n.edge,1)
            if !n.edge[i].hybrid push!(hyb,n.edge[i]); end
        end
        #println("hyb tiene $([e.number for e in hyb])")
        for e in hyb
            #println("se va a borrar a $(e.number)")
            removeEdge!(n,e);
            removeNode!(n,e);
            setEdge!(n1,e);
            setNode!(e,n1);
        end
        #println("now node $(n1.number) has the edges $([e.number for e in n1.edge])")
        setEdge!(n,ed1);
        setNode!(ed1,[n,n1]);
        pushNode!(net,n1);
        pushEdge!(net,ed1);
        if size(n1.edge,1) > 3
            solvePolytomy!(net,n1);
        end
    else
        error("cannot expand children of a tree node.")
    end
end

# function to clean topology after readTopology
# looks for:
# TREE:
# - all tree edges must have gamma=1. fixit: cannot point out which doesn't,
#   only shows error.
# - internal nodes with only 2 edges and solves this case.
# - polytomies and choose one resolution at random, issuing a warning
# NETWORK:
# - number of hybrid edges per hybrid node:
#   if 0,1: error (with warning in old functions)
#   if >2: error of hybrid polytomy
#   if 2: check number of tree edges
# - number of tree edges per hybrid node:
#   if 0: leaf hybrid, add child
#   if >1: expand child
#   if 1: check values of gamma:
# - gammas: need to sum to one and be present.
#   error if they do not sum up to one
#   default values of 0.1,0.9 if not present
# leaveRoot=true: leaves the root even if it has only 2 edges (for plotting), default=false
function cleanAfterRead!(net::HybridNetwork, leaveRoot::Bool)
    # set e.containRoot to !e.hybrid: updated later by updateAllReadTopology as required by snaq!
    for e in net.edge e.containRoot = !e.hybrid; end
    nodes = copy(net.node)
    for n in nodes
        if isNodeNumIn(n,net.node) # very important to check
            if size(n.edge,1) == 2
                if !n.hybrid
                    if !leaveRoot || !isEqual(net.node[net.root],n) #if n is the root
                        deleteIntNode!(net,n);
                    end
                else
                    hyb = count([e.hybrid for e in n.edge]);
                    if hyb == 1
                        deleteIntNode!(net,n);
                    end
                end
            end
            if !n.hybrid
                if size(n.edge,1) > 3
                    @debug "warning: polytomy found in node $(n.number), random resolution chosen"
                    solvePolytomy!(net,n);
                end
                hyb = count([e.hybrid for e in n.edge])
                if hyb == 1
                    n.hasHybEdge == true;
                elseif hyb > 1
                    @warn "strange tree node $(n.number) with more than one hybrid edge, intersecting cycles maybe"
                end
            else
                hyb = count([e.hybrid for e in n.edge]);
                tre = length(n.edge) - hyb
                if hyb > 2
                    error("hybrid node $(n.number) has more than two hybrid edges attached to it: polytomy that cannot be resolved without intersecting cycles.")
                elseif hyb == 1
                    hybnodes = count([n.hybrid for n in net.node]);
                    if hybnodes == 1
                        error("only one hybrid node number $(n.number) with name $(net.names[n.number]) found with one hybrid edge attached")
                    else
                        error("current hybrid node $(n.number) with name $(net.names[n.number]) has only one hybrid edge attached. there are other $(hybnodes-1) hybrids out there but this one remained unmatched")
                    end
                elseif hyb == 0
                    @warn "hybrid node $(n.number) is not connected to any hybrid edges, it was transformed to tree node"
                    n.hybrid = false;
                else # 2 hybrid edges
                    if tre == 0 #hybrid leaf
                        @warn "hybrid node $(n.number) is a leaf, so we add an extra child"
                        addChild!(net,n);
                    elseif tre > 1
                        @warn "hybrid node $(n.number) has more than one child so we need to expand with another node"
                        expandChild!(net,n);
                    end
                    suma = sum([e.hybrid ? e.gamma : 0.0 for e in n.edge]);
                    # synchronizePartnersData! already made suma ≈ 1.0, when non-missing,
                    # and already synchronized isMajor, even when γ's ≈ 0.5
                    if suma == -2.0 # hybrid edges have no gammas in newick description
                        println("hybrid edges for hybrid node $(n.number) have missing gamma's, set default: 0.9,0.1")
                        for e in n.edge
                            if e.hybrid
                                e.gamma = (e.isMajor ? 0.9 : 0.1)
                            end
                        end
                    end
                end
            end
        end
    end
    for e in net.edge
        if e.hybrid
          if e.gamma < 0.0 || e.gamma > 1.0 # no -1.0 (missing) γ's by now
            error("hybrid edge $(e.number) with γ = $(e.gamma) outside of [0,1]")
          end
        else
          e.gamma == 1.0 ||
            error("tree edge $(e.number) with γ = $(e.gamma) instead of 1.0")
        end
    end
end

cleanAfterRead!(net::HybridNetwork) = cleanAfterRead!(net,false)

# function to search for the hybrid nodes in a read network after cleaning it
# and store this information as a network's attribute
function storeHybrids!(net::HybridNetwork)
    flag = true;
    hybrid = nothing
    try
        hybrid = searchHybridNode(net)
    catch
        #@warn "topology read is a tree as it has no hybrid nodes"
        flag = false;
    end
    if(flag)
        net.hybrid = hybrid;
        net.numHybrids = size(hybrid,1);
    end
    return nothing
end

# function to update the read topology after reading
# it will go over the net.hybrid array and check each
# of the hybridization events defined to update:
# - in cycle
# - contain root
# - gammaz
# it uses updateAllNewHybrid! function that
# returns: success (bool), hybrid, flag, nocycle, flag2, flag3
# if tree read, also check that contain root is true for all, ishybrid and hashybedge is false
# warning: needs to have run storeHybrids! before
# warning: it will stop when finding one conflicting hybrid
function updateAllReadTopology!(net::HybridNetwork)
    if(isTree(net))
        #@warn "not a network read, but a tree as it does not have hybrid nodes"
        all((e->e.containRoot), net.edge) ? nothing : error("some tree edge has contain root as false")
        all((e->!e.hybrid), net.edge) ? nothing : error("some edge is hybrid and should be all tree edges in a tree")
        all((n->!n.hasHybEdge), net.node) ? nothing : error("some tree node has hybrid edge true, but it is a tree, there are no hybrid edges")
    else
        if(!net.cleaned)
            for n in net.hybrid
                success,hyb,flag,nocycle,flag2,flag3 = updateAllNewHybrid!(n,net,false,true,false)
                if(!success)
                    @warn "current hybrid $(n.number) conflicts with previous hybrid by intersecting cycles: $(!flag), nonidentifiable topology: $(!flag2), empty space for contain root: $(!flag3), or does not create a cycle (probably problem with the root placement): $(nocycle)."
                    #net.cleaned = false
                end
            end
            @debug "before update partition"
            @debug begin printPartitions(net); "printed partitions" end
            for n in net.hybrid #need to updatePartition after all inCycle
                nocycle, edgesInCycle, nodesInCycle = identifyInCycle(net,n);
                updatePartition!(net,nodesInCycle)
                @debug begin printPartitions(net);
                    "partitions after updating partition for hybrid node $(n.number)"
                end
            end
        end
    end
end

# cleanAfterReadAll includes all the step to clean a network after read
function cleanAfterReadAll!(net::HybridNetwork, leaveRoot::Bool)
    @debug "check for 2 hybrid edges at each hybrid node -----"
    check2HybEdges(net)
    @debug "cleanBL -----"
    cleanBL!(net)
    @debug "cleanAfterRead -----"
    cleanAfterRead!(net,leaveRoot)
    @debug "updateAllReadTopology -----"
    updateAllReadTopology!(net) #fixit: it could break if leaveRoot = true (have not checked it), but we need to updateContainRoot
    if(!leaveRoot)
        @debug "parameters -----"
        parameters!(net)
    end
    @debug "check root placement -----"
    checkRootPlace!(net)
    net.node[net.root].leaf && @warn "root node $(net.node[net.root].number) is a leaf, so when plotting net, it can look weird"
    net.cleaned = true #fixit: set to false inside updateAllReadTopology if problem encountered
    net.isRooted = false
end

cleanAfterReadAll!(net::HybridNetwork) = cleanAfterReadAll!(net,false)

# function to read a topology from file name/tree directly and update it
# by calling updateAllReadTopology after
# leaveRoot=true if the root will not be deleted even if it has only 2 edges
# used for plotting (default=false)
# warning: if leaveRoot=true, net should not be used outside plotting, things will crash
function readTopologyUpdate(file::AbstractString, leaveRoot::Bool,verbose::Bool)
    @debug "readTopology -----"
    net = readTopology(file,verbose)
    cleanAfterReadAll!(net,leaveRoot)
    return net
end

readTopologyUpdate(file::AbstractString) = readTopologyUpdate(file, false, true)
readTopologyUpdate(file::AbstractString,verbose::Bool) = readTopologyUpdate(file, false, verbose)

"""
    readTopologyLevel1(filename)
    readTopologyLevel1(parenthetical format)

same as readTopology, reads a tree or network from parenthetical
format, but this function enforces the necessary conditions for any
starting topology in SNaQ: non-intersecting cycles, no polytomies,
unrooted. It sets any missing branch length to 1.0.

If the network has a bad diamond II (in which edge lengths are γ's are not identifiable)
and if the edge below this diamond has a length `t` different from 0, then this length is
set back to 0 and the major parent hybrid edge is lengthened by `t`.
"""
readTopologyLevel1(file::AbstractString) = readTopologyUpdate(file, false, true)

# aux function to check if the root is placed correctly, and re root if not
# warning: it needs updateContainRoot set
function checkRootPlace!(net::HybridNetwork; verbose=false::Bool, outgroup="none"::AbstractString)
    if(outgroup == "none")
        if(!canBeRoot(net.node[net.root]))
            verbose && println("root node $(net.node[net.root].number) placement is not ok, we will change it to the first found node that agrees with the direction of the hybrid edges")
            for i in 1:length(net.node)
                if(canBeRoot(net.node[i]))
                    net.root = i
                    break
                end
            end
        end
    else # outgroup
        tmp = findall(n -> n.name == outgroup, net.leaf)
        if length(tmp)==0
            error("leaf named $(outgroup) was not found in the network.")
        elseif length(tmp)>1
            error("several leaves were found with name $(outgroup).")
        end
        leaf = net.leaf[tmp[1]]
        leaf.leaf || error("found outgroup not a leaf: $(leaf.number), $(outgroup)")
        length(leaf.edge) == 1 || error("found leaf with more than 1 edge: $(leaf.number)")
        other = getOtherNode(leaf.edge[1],leaf);
        if(canBeRoot(other))
            net.root = getIndexNode(other.number,net)
        else
            throw(RootMismatch("outgroup $(outgroup) contradicts direction of hybrid edges"))
        end
    end
    canBeRoot(net.node[net.root]) || error("tried to place root, but couldn't. root is node $(net.node[net.root])")
end

"""
    writeSubTree!(IO, network, dendroscope::Bool, namelabel::Bool,
                  round_branch_lengths::Bool, digits::Integer,
                  internallabel::Bool)

Write to IO the extended newick format (parenthetical description)
of a network.
If written for dendroscope, inheritance γ's are not written.
If `namelabel` is true, taxa are labelled by their names;
otherwise taxa are labelled by their number IDs.
If unspecified, branch lengths and γ's are rounded to 3 digits.
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)
    rootnode = net.node[net.root]
    if net.numNodes > 1
        print(s,"(")
        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

"""
    writeSubTree!(IO, node, edge, dendroscope::Bool, namelabel::Bool,
                  round_branch_lengths::Bool, digits::Integer, internallabel::Bool)

Write the extended newick format of the sub-network rooted at
`node` and assuming that `edge` is a parent of `node`.

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:

```julia
net = readTopology("(((A,(B)#H1:::0.9),(C,#H1:::0.1)),D);")
directEdges!(net)
s = IOBuffer()
writeSubTree!(s, net.node[7], nothing, false, true)
String(take!(s))
```

Used by [`writeTopology`](@ref).
"""
writeSubTree!(s,n,parent,di,namelabel) =
    writeSubTree!(s,n,parent,di,namelabel, true,3,true)

# "parent' is assumed to be adjancent to "node". not checked.
# algorithm comes from "parent": do not traverse again.
function writeSubTree!(s::IO, n::Node, parent::Union{Edge,Nothing},
    di::Bool, namelabel::Bool, roundBL::Bool, digits::Integer, internallabel::Bool)
    # subtree below node n:
    if !n.leaf && (parent == nothing || parent.isMajor) # do not descent below a minor hybrid edge
        print(s,"(")
        firstchild = true
        for e in n.edge
            e != parent || continue # skip parent edge where we come from
            if parent == nothing    # skip if n = child of e
                n != getChild(e) || continue
            end
            (e.hybrid && getChild(e)==n) && continue # no going up minor hybrid
            firstchild || print(s, ",")
            firstchild = false
            child = getOtherNode(e,n)
            writeSubTree!(s,child,e, di,namelabel, roundBL, digits, internallabel)
        end
        print(s,")")
    end
    # node label:
    if parent != nothing && parent.hybrid
        print(s, "#")
        print(s, (namelabel ? n.name : string("H", n.number)))
        n.name != "" || parent.isMajor || @warn "hybrid node $(n.number) has no name"
    elseif internallabel || n.leaf
        print(s, (namelabel ? n.name : n.number))
    end
    # branch lengths and γ, if available:
    printBL = false
    if parent != nothing && parent.length != -1.0 # -1.0 means missing
        print(s,string(":",(roundBL ? round(parent.length, digits=digits) : parent.length)))
        printBL = true
    end
    if parent != nothing && parent.hybrid && !di # && (!printID || !n.isBadDiamondI))
        if(parent.gamma != -1.0)
            if(!printBL) print(s,":"); end
            print(s,string("::",(roundBL ? round(parent.gamma, digits=digits) : parent.gamma)))
        end
    end
    if parent == nothing
        print(s, ";")
    end
end

# see full docstring below
# Need HybridNetwork input, since QuartetNetwork does not have root.
function writeTopologyLevel1(net0::HybridNetwork, di::Bool, str::Bool, namelabel::Bool,outgroup::AbstractString, printID::Bool, roundBL::Bool, digits::Integer, multall::Bool)
    s = IOBuffer()
    writeTopologyLevel1(net0,s,di,namelabel,outgroup,printID,roundBL,digits, multall)
    if str
        return String(take!(s))
    else
        return s
    end
end

# warning: I do not want writeTopologyLevel1 to modify the network if outgroup is given! thus, we have updateRoot, and undoRoot
function writeTopologyLevel1(net0::HybridNetwork, s::IO, di::Bool, namelabel::Bool,
           outgroup::AbstractString, printID::Bool, roundBL::Bool, digits::Integer, multall::Bool)
    global CHECKNET
    net = deepcopy(net0) #writeTopologyLevel1 needs containRoot, but should not alter net0
    # net.numBad == 0 || println("network with $(net.numBad) bad diamond I. Some γ and edge lengths t are not identifiable, although their γ * (1-exp(-t)) are.")
    if printID
        setNonIdBL!(net) # changes non identifiable BL to -1.0, except those in/below bad diamonds/triangles.
    end
    assignhybridnames!(net)
    if(net.numNodes == 1)
        print(s,string(net.node[net.root].number,";")) # error if 'string' is an argument name.
    else
        if(!isTree(net) && !net.cleaned)
            @debug "net not cleaned inside writeTopologyLevel1, need to run updateContainRoot"
            for n in net.hybrid
                flag,edges = updateContainRoot!(net,n)
                flag || error("hybrid node $(n.hybrid) has conflicting containRoot")
            end
        end
        updateRoot!(net,outgroup)
        #@debug begin printEverything(net); "printed everything" end
        CHECKNET && canBeRoot(net.node[net.root])
        if(multall)
            mergeLeaves!(net)
        end
        writeSubTree!(s, net, di,namelabel, roundBL,digits,false)
    end
    # outgroup != "none" && undoRoot!(net) # not needed because net is deepcopy of net0
    # to delete 2-degree node, for snaq.
end

writeTopologyLevel1(net::HybridNetwork,di::Bool,str::Bool,namelabel::Bool,outgroup::AbstractString,printID::Bool) = writeTopologyLevel1(net,di,str,namelabel,outgroup,printID, false,3, false)
# above: default roundBL=false (at unused digits=3 decimal places)
writeTopologyLevel1(net::HybridNetwork,printID::Bool) = writeTopologyLevel1(net,false, true,true,"none",printID, false, 3, false)
writeTopologyLevel1(net::HybridNetwork,outgroup::AbstractString) = writeTopologyLevel1(net,false, true,true,outgroup,true, false, 3, false)
writeTopologyLevel1(net::HybridNetwork,di::Bool,outgroup::AbstractString) = writeTopologyLevel1(net,di, true,true,outgroup,true, false, 3, false)

"""
`writeTopologyLevel1(net::HybridNetwork)`

Write the extended Newick parenthetical format of a
level-1 network object with many optional arguments (see below).
Makes a deep copy of net: does *not* modify `net`.

- di=true: write in format for Dendroscope (default false)
- namelabel=true: If `namelabel` is true, taxa are labelled by their names;
otherwise taxa are labelled by their numbers (unique identifiers).
- outgroup (string): name of outgroup to root the tree/network.
  if "none" is given, the root is placed wherever possible.
- printID=true, only print branch lengths for identifiable egdes
  according to the snaq estimation procedure (default false)
  (true inside of `snaq!`.)
- round: rounds branch lengths and heritabilities γ (default: true)
- digits: digits after the decimal place for rounding (defult: 3)
- string: if true (default), returns a string,
  otherwise returns an IOBuffer object.
- multall: (default false). set to true when there are multiple
  alleles per population.

The topology may be written using a root different than net.root,
if net.root is incompatible with one of more hybrid node.
Missing hybrid names are written as "#Hi" where "i" is the hybrid node number if possible.
""" #"
writeTopologyLevel1(net::HybridNetwork; di=false::Bool, string=true::Bool, namelabel=true::Bool,
    outgroup="none"::AbstractString, printID=false::Bool, round=false::Bool, digits=3::Integer,
    multall=false::Bool) =
writeTopologyLevel1(net, di, string, namelabel, outgroup, printID, round, digits, multall)

# function to check if root is well-placed
# and look for a better place if not
# searches on net.node because net.root is the index in net.node
# if we search in net.edge, we then need to search in net.node
function updateRoot!(net::HybridNetwork, outgroup::AbstractString)
    checkroot = false
    if(outgroup == "none")
        @debug "no outgroup defined"
        checkroot = true
    else
        println("outgroup defined $(outgroup)")
        index = findfirst(n -> outgroup == n.name, net.node)
        index != nothing ||
            error("outgroup $(outgroup) not in net.names $(net.names)")
        node = net.node[index]
        node.leaf || error("outgroup $(outgroup) is not a leaf in net")
        length(net.node[index].edge) == 1 || error("strange leaf $(outgroup), node number $(net.node[index].number) with $(length(net.node[index].edge)) edges instead of 1")
        edge = net.node[index].edge[1]
        if(edge.containRoot)
            DEBUGC && @debug "creating new node in the middle of the external edge $(edge.number) leading to outgroup $(node.number)"
            othernode = getOtherNode(edge,node)
            removeEdge!(othernode,edge)
            removeNode!(othernode,edge)
            max_edge = maximum([e.number for e in net.edge]);
            max_node = maximum([e.number for e in net.node]);
            newedge = Edge(max_edge+1) #fixit: maybe this edge not identifiable, need to add that check
            newnode = Node(max_node+1,false,false,[edge,newedge])
            if(net.cleaned && !isTree(net) && !isempty(net.partition)) # fixit: this will crash if network estimated with snaq, and then manipulated
                part = whichPartition(net,edge)
                push!(net.partition[part].edges,newedge)
            end
            setNode!(edge,newnode)
            setNode!(newedge,newnode)
            setEdge!(othernode,newedge)
            setNode!(newedge,othernode)
            pushEdge!(net,newedge)
            pushNode!(net,newnode)
            t = edge.length
                setLength!(edge,t/2)
                setLength!(newedge,t/2)
            net.root = length(net.node) #last node is root
       else
            @warn "external edge $(net.node[index].edge[1].number) leading to outgroup $(outgroup) cannot contain root, root placed wherever"
            checkroot = true
       end
    end
    if(checkroot && !isTree(net))
        checkRootPlace!(net)
    end
 end

# function to check if a node could be root
# by the containRoot attribute of edges around it
function canBeRoot(n::Node)
    !n.hybrid || return false
    #!n.hasHybEdge || return false #need to allow for some reason, check ipad notes
    !n.leaf || return false
    return any([e.containRoot for e in n.edge])
end

# function to delete the extra node created in updateRoot
# this extra node is needed to be able to compare networks with the distance function
# but if left in the network, everything crashes (as everything assumes three edges per node)
# fromUpdateRoot=true if called after updateRoot (in which case leaf has to be a leaf), ow it is used in readTopUpdate
function undoRoot!(net::HybridNetwork, fromUpdateRoot::Bool)
    if(length(net.node[net.root].edge) == 2)
        root = net.node[net.root]
        leaf = getOtherNode(root.edge[1],root).leaf ? getOtherNode(root.edge[1],root) : getOtherNode(root.edge[2],root)
        (fromUpdateRoot && leaf.leaf) || error("root should have one neighbor leaf which has to be the outgroup defined")
        deleteIntLeafWhile!(net,root,leaf);
    end
end

undoRoot!(net::HybridNetwork) = undoRoot!(net, true)

# .out file from snaq written by optTopRuns
"""
    readSnaqNetwork(output file)

Read the estimated network from a `.out` file generated by `snaq!`.
The network score is read also, and stored in the network's field `.loglik`.

Warning: despite the name "loglik", this score is only proportional
to the network's pseudo-deviance: the lower, the better.
Do NOT use this score to calculate an AIC or BIC (etc.) value.
"""
function readSnaqNetwork(file::AbstractString)
    open(file) do s
        line = readline(s)
        line[1] == '(' ||
          error("output file $(file) does not contain a tree in the first line, instead it has $(line); or we had trouble finding ploglik.")
        # println("Estimated network from file $(file): $(line)")
        net = readTopology(line)
        # readTopologyUpdate is inadequate: would replace missing branch lengths, which are unidentifiable, by 1.0 values
        try
            vec = split(line,"-Ploglik = ")
            net.loglik = parse(Float64,vec[2])
        catch e
            @warn "could not find the network score; the error was:"
            rethrow(e)
        end
        return net
    end
end

# function to change negative branch lengths to 1.0 for starting topology
# and to change big branch lengths to 10.0
# also uses setLength for all edges
function cleanBL!(net::HybridNetwork)
    ##println("missing branch lengths will be set to 1.0")
    for e in net.edge
        if(e.length < 0)
            setLength!(e,1.0)
        elseif(e.length > 10.0)
            setLength!(e,10.0)
        else
            setLength!(e,e.length)
        end
    end
end

# function to read multiple topologies
# - calls readInputTrees in readData.jl, which
#   calls readTopologyUpdate here, for level 1 networks.
# - read a file and create one object per line read
# (each line starting with "(" will be considered a topology)
# the file can have extra lines that are ignored
# returns an array of HybridNetwork objects (that can be trees)
function readMultiTopologyLevel1(file::AbstractString)
    readInputTrees(file)
end

"""
`readMultiTopology(file)`

Read a text file with a list of networks in parenthetical format (one per line).
Each network is read with [`readTopology`](@ref).
Return an array of HybridNetwork object.
"""
function readMultiTopology(file::AbstractString)
    s = open(file)
    numl = 1
    vnet = HybridNetwork[];
    for line in eachline(s)
        line = strip(line) # remove spaces
        c = isempty(line) ? "" : line[1]
        if(c == '(')
           try
               push!(vnet, readTopology(line,false)) # false for non-verbose
           catch err
               print("skipped phylogeny on line $(numl) of file $file: ")
               if :msg in fieldnames(typeof(err)) println(err.msg); else println(typeof(err)); end
           end
        end
        numl += 1
    end
    close(s)
    return vnet
end

"""
    writeMultiTopology(nets, file_name; append=false)
    writeMultiTopology(nets, IO)

Write an array of networks in parenthetical extended Newick format, one network per line.
Use the option append=true to append to the file. Otherwise, the default is to create a new
file or overwrite it, if it already existed.
Each network is written with `writeTopology`.

# Examples
```
julia> net = [readTopology("(D,((A,(B)#H7:::0.864):2.069,(F,E):3.423):0.265,(C,#H7:::0.1361111):10);"),
              readTopology("(A,(B,C));"),readTopology("(E,F);"),readTopology("(G,H,F);")];

julia> writeMultiTopology(net, "fournets.net") # to (over)write to file "fournets.net"
julia> writeMultiTopology(net, "fournets.net", append=true) # to append to this file
julia> writeMultiTopology(net, stdout)         # to write to the screen (standard out)
(D,((A,(B)#H7:::0.864):2.069,(F,E):3.423):0.265,(C,#H7:::0.1361111):10.0);
(A,(B,C));
(E,F);
(G,H,F);
```
"""
function writeMultiTopology(n::Vector{HybridNetwork},file::AbstractString; append::Bool=false)
    mode = (append ? "a" : "w")
    open(file, mode) do s
    writeMultiTopology(n,s)
    end # closes file safely
end

function writeMultiTopology(net::Vector{HybridNetwork},s::IO)
    for i in 1:length(net)
      try
        # writeTopologyLevel1(net[i],s,false,true,"none",false,false,3)
        writeTopology(net[i],s) # no rounding, not for dendroscope
        write(s,"\n")
      catch err
        if isa(err, RootMismatch) # continue writing other networks in list
            @error "\nError with topology $i:\n" * err.msg
        else rethrow(err); end
      end
    end
end

"""
    writeTopology(net)
    writeTopology(net, filename)
    writeTopology(net, IO)

Write the parenthetical extended Newick format of a network,
as a string, to a file or to an IO buffer / stream.
Optional arguments (default values):

- di (false): write in format for Dendroscope
- round (false): rounds branch lengths and heritabilities γ
- digits (3): digits after the decimal place for rounding
- append (false): if true, appends to the file
- internallabel (true): if true, writes internal node labels

If the current root placement is not admissible, other placements are tried.
The network is updated with this new root placement, if successful.

Uses lower-level function [`writeSubTree!`](@ref).
"""
function writeTopology(n::HybridNetwork, file::AbstractString; append::Bool=false,
        round=false::Bool, digits=3::Integer, di=false::Bool, # keyword arguments
        internallabel=true::Bool)
    mode = (append ? "a" : "w")
    s = open(file, mode)
    writeTopology(n,s,round,digits,di,internallabel)
    write(s,"\n")
    close(s)
end

function writeTopology(n::HybridNetwork;
        round=false::Bool, digits=3::Integer, di=false::Bool, # keyword arguments
        internallabel=true::Bool)
    s = IOBuffer()
    writeTopology(n,s,round,digits,di,internallabel)
    return String(take!(s))
end

function writeTopology(net::HybridNetwork, s::IO,
        round=false::Bool, digits=3::Integer, di=false::Bool, # optional arguments
        internallabel=true::Bool)
    # check/find admissible root: otherwise could be trapped in infinite loop
    rootsaved = net.root
    changeroot = false
    msg = ""
    try
        directEdges!(net)
    catch err
        if isa(err, RootMismatch)
            println(err.msg * "\nCannot write topology with current root.")
            changeroot = true
        else rethrow(err); end
    end
    while changeroot
        for e in net.edge
          # parents of hybrid edges should be sufficient, but gives weird look
          #if e.hybrid
            i = getIndex(getParent(e), net)
            net.root = i
            try
                directEdges!(net)
                print("Setting root at node $(net.node[i].number) (net.root = $i)\n\n")
                print(msg)
                changeroot = false
                break # stop loop over edges
            catch err
                if !isa(err, RootMismatch) rethrow(err); end
            end
          #end
        end
        if changeroot # none of hybrid edges worked
            net.root = rootsaved
            throw(RootMismatch("Could not find admissible root. Cannot write topology."))
            changeroot=false # safety exit of while (but useless)
        end
    end
    # finally, write parenthetical format
    writeSubTree!(s,net,di,true,round,digits,internallabel)
    # namelabel = true: to print leaf & node names (labels), not numbers
end

"""
    hybridlambdaformat(net::HybridNetwork; prefix="I")

Output `net` as a string in the format that the
[Hybrid-Lambda](https://github.com/hybridLambda/hybrid-Lambda)
simulator expects, namely:
- all internal nodes are named, including the root
- hybrid nodes are written as `H6#γ1:length1` and `H6#γ1:length2`
  instead of `#H6:length1::γ1` and `#H6:length2::γ2`
  (note the samme γ value expected by Hybrid-Lambda)

This is a modified version of the
[extended Newick](https://doi.org/10.1186/1471-2105-9-532) format.

Optional keyword argument: `prefix`, "I" by default.
Internal nodes without a name will be given names like "I1", "I2", etc.
See [`nameinternalnodes!`](@ref) to add node names (modifies the input network).

# examples

```jldoctest
julia> net = readTopology("((a:1,(b:1)#H1:1::0.8):5,(#H1:0::0.2,c:1):1);");

julia> hybridlambdaformat(net) # net is unchanged here
"((a:1.0,(b:1.0)H1#0.8:1.0)I1:5.0,(H1#0.8:0.0,c:1.0)I2:1.0)I3;"

julia> # using PhyloPlots; plot(net, :R, showNodeNumber=true) # shows that node -2 is the root

julia> rotate!(net, -2)

julia> writeTopology(net) # now the minor edge with γ=0.2 appears first
"((#H1:0.0::0.2,c:1.0):1.0,(a:1.0,(b:1.0)#H1:1.0::0.8):5.0);"

julia> hybridlambdaformat(net)
"((H1#0.2:0.0,c:1.0)I2:1.0,(a:1.0,(b:1.0)H1#0.2:1.0)I1:5.0)I3;"

julia> net = readTopology("((((B)#H1:::.6)#H2,((D,C,#H2:::0.8),(#H1,A))));"); # 2 reticulations, no branch lengths

julia> writeTopology(net, round=true)
"(#H2:::0.2,((D,C,((B)#H1:::0.6)#H2:::0.8),(#H1:::0.4,A)));"

julia> hybridlambdaformat(net; prefix="int")
"(H2#0.2,((D,C,((B)H1#0.6)H2#0.2)int1,(H1#0.6,A)int2)int3)int4;"
```
"""
function hybridlambdaformat(net::HybridNetwork; prefix="I")
  net = deepcopy(net) # binding to new object
  nameinternalnodes!(net, prefix)
  str1 = writeTopology(net, round=true, digits=15) # internallabels=true by default
  rx_noBL = r"#(H[\w\d]+)::\d*\.?\d*:(\d*\.?\d*)"
  subst_noBL = s"\1#\2"
  rx_withBL = r"#(H[\w\d]+):(\d*\.?\d*):\d*\.?\d*:(\d*\.?\d*)"
  subst_withBL = s"\1#\3:\2"
  str2 = replace(replace(str1, rx_noBL => subst_noBL), rx_withBL => subst_withBL)
  ## next: replace the γ2 of the second occurrence by γ1 from the first occurrence:
  ## this is what Hybrid-Lambda wants...
  nh = length(net.hybrid)
  m = eachmatch(r"H[^#]*#", str2)
  hboth = collect(h.match for h in m)
  hone = unique(hboth)
  length(hboth) == 2length(hone) || error("did not find Hname# twice for some one (or more) of the hybrid names.")
  str3 = str2
  for hname in hone
    rx = Regex( hname * raw"(?<gamma1>\d*\.?\d*)(?<middle>.*)" * hname * raw"(?<gamma2>\d*\.?\d*)")
    subst = SubstitutionString(hname * raw"\g<gamma1>\g<middle>" * hname * raw"\g<gamma1>")
    str3 = replace(str3, rx => subst)
  end
  return str3
end

"""
    nameinternalnodes!(net::HybridNetwork, prefix)

Add names to nodes in `net` that don't already have a name.
Leaves already have names; but if not, they will be give names as well.
New node names will be of the form "prefixI" where I is an integer.

# examples
```jldoctest
julia> net = readTopology("((a:1,(b:1)#H1:1::0.8):5,(#H1:0::0.2,c:1):1);");

julia> PhyloNetworks.nameinternalnodes!(net, "I") # by default, shown without internal node names
HybridNetwork, Rooted Network
7 edges
7 nodes: 3 tips, 1 hybrid nodes, 3 internal tree nodes.
tip labels: a, b, c
((a:1.0,(b:1.0)#H1:1.0::0.8):5.0,(#H1:0.0::0.2,c:1.0):1.0);

julia> writeTopology(net) # by default, writeTopology shows internal names if they exist
"((a:1.0,(b:1.0)#H1:1.0::0.8)I1:5.0,(#H1:0.0::0.2,c:1.0)I2:1.0)I3;"

julia> net = readTopology("((int5:1,(b:1)#H1:1::0.8):5,(#H1:0::0.2,c:1):1);"); # one taxon name starts with "int"

julia> PhyloNetworks.nameinternalnodes!(net, "int");

julia> writeTopology(net)
"((int5:1.0,(b:1.0)#H1:1.0::0.8)int6:5.0,(#H1:0.0::0.2,c:1.0)int7:1.0)int8;"
```
"""
function nameinternalnodes!(net::HybridNetwork, prefix)
  # get maximum index I of nodes whose names are already like: prefixI
  rx = Regex("^$(prefix)(\\d+)\$")
  nexti = 1
  for node in net.node
    node.name != "" || continue # skip nodes with empty names
    m = match(rx, node.name)
    m != nothing || continue
    nexti = max(nexti, parse(Int, m.captures[1])+1)
  end
  # assign names like: prefixI for I = nexti, nexti+1, etc.
  for node in net.node
    node.name == "" || continue # skip nodes with non-empty names
    node.name = prefix * string(nexti)
    nexti += 1
  end
  return net
end

###############################################################################
## Generate symetric tree
###############################################################################

"""
    symmetricTree(n, i=1)

Create a string with a symmetric tree with 2^n tips, numbered from i to i+2^n-1.
All the branch length are set equal to 1.
The tree can be created with function readTopology.
"""
function symmetricTree(n::Int, ell::Real, i=1::Int)
    # Build tree
    tree = "A$(i-1+2^n):$(ell)"
    if n==0 return("("*"A$(i-1+2^n):0"*");") end
    for k in 1:(n-1)
        tree = "(" * tree * "," * tree * "):$(ell)"
    end
    tree = "(" * tree * "," * tree * ");"
    # Rename tips
    for k in (2^n-1):-1:1
        tree = replace(tree, "A$(i-1+k+1):$(ell)"=>"A$(i-1+k):$(ell)", count=k)
    end
    return(tree)
end

"""
    symmetricNet(n, i, j, gamma)

Create a string with a symmetric net with 2^n tips, numbered from 1 to 2^n
All the branch length are set equal to 1.
One hybrid branch, going from level i to level j is added, with weigth gamma.
The tree can be created with function readTopology.
"""
function symmetricNet(n::Int, i::Int, j::Int, gamma::Real, ell::Real)
    # Checks
    if (n < i || i < j || j < 1) error("Must be n > i > j > 0") end
    # Underlying tree
    tree = symmetricTree(n, ell)
    ## start hyb
    op = "(A1:$(ell)"
    clo = "A$(2^(i-1)):$(ell)):$(ell)"
    for k in 3:i
        op = "("*op
        clo = clo*"):$(ell)"
    end
    clobis = clo[1:(length(clo)-length("$(ell)"))]*"$(0.5*ell)"
    tree = replace(tree,  op => "(#H:$(ell)::$(gamma),"*op)
    tree = replace(tree, clo => clobis*"):$(0.5*ell)")
    ## end hyb
    op = "A$(2^(i-1)+1):$(ell)"
    clo = "A$(2^(i-1) + 2^(j-1)):$(ell)"
    for k in 2:j
        op = "("*op
        clo = clo*"):$(ell)"
    end
    clobis = clo[1:(length(clo)-length("$(ell)"))]*"$(0.5*ell)"
    tree = replace(tree, op  => "("*op)
    tree = replace(tree, clo => clobis*")#H:$(0.5*ell)::$(1-gamma)")
    return(tree)
end

"""
    symmetricNet(n, h, gamma)

Create a string with a symmetric net with 2^n tips, numbered from 1 to 2^n
The total height of the network is set to 1.
Hybrids are added from level h to h-1 symmetrically.
"""
function symmetricNet(n::Int, h::Int, gamma::Real, i=1::Int)
    # Checks
    if (n < h || h < 2) error("Must be n >= h > 1.") end
    # length of branch
    ell = 1.0/n
    # Element net
    net = symmetricNet(h, h, h-1, gamma, ell)
    # Iterate
    if n==h return(net) end
    net = chop(net)
    for k in 1:(n-h)
        net = "(" * net * ":$(ell)," * net * ":$(ell))"
    end
    net = net * ";"
    # Rename hybrids
    net = replace(net, "H" => "H$(2^(n-h))")
    for k in (2^(n-h)-1):-1:1
        net = replace(net, "H$(k+1)" => "H$(k)", count=2*k)
    end
    # Rename tips
    for k in (2^h):-1:1
        net = replace(net, "A$(k):" => "A$(i-1+2^n):")
    end
    for k in (2^n-1):-1:1
        net = replace(net, "A$(i-1+k+1):" => "A$(i-1+k):", count=k)
    end
    return(net)
end


text:
# ------------------------------------------------------------------
# Copyright (c) 2017, Júlio Hoffimann Mendes <juliohm@stanford.edu>
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

function digest(solution::EstimationSolution{<:RegularGrid})
  # solution variables
  variables = keys(solution.mean)

  # get the size of the grid
  sz = size(solution.domain)

  # build dictionary pairs
  pairs = []
  for var in variables
    M = reshape(solution.mean[var], sz)
    V = reshape(solution.variance[var], sz)

    push!(pairs, var => Dict(:mean => M, :variance => V))
  end

  # output dictionary
  Dict(pairs)
end

ext: jl


text:
module Reinforce

using Reexport
@reexport using StatsBase
using Distributions
using RecipesBase

@reexport using LearnBase
import LearnBase: learn!, transform!, grad!, params, grad

using LearningStrategies
import LearningStrategies: setup!, hook, finished, cleanup!

export
  AbstractEnvironment,
  reset!,
  step!,
  reward,
  state,
  finished,
  actions,
  ismdp,
  maxsteps,

  AbstractPolicy,
  RandomPolicy,
  OnlineGAE,
  OnlineActorCritic,
  EpisodicActorCritic,
  action,

  AbstractState,
  StateVector,
  History,
  state!,

  Episode,
  Episodes,
  run_episode

# ----------------------------------------------------------------
# Implement this interface for a new environment

abstract type AbstractEnvironment end

"""
    reset!(env) -> env

Reset an environment.
"""
function reset! end

"""
    r, s′ = step!(env, s, a)

Move the simulation forward, collecting a reward and getting the next state.
"""
function step! end

# note for developers: you should also implement Base.done(env) for episodic environments
finished(env::AbstractEnvironment, s′) = false

"""
    A = actions(env, s)

Return a list/set/description of valid actions from state `s`.
"""
function actions end

# note for developers: you don't need to implement these if you have state/reward fields

"""
    s = state(env)

Return the current state of the environment.
"""
state(env::AbstractEnvironment) = env.state

"""
    r = reward(env)

Return the current reward of the environment.
"""
reward(env::AbstractEnvironment) = env.reward

"""
    ismdp(env)::Bool

An environment may be fully observable (MDP) or partially observable (POMDP).
In the case of a partially observable environment,
the state `s` is really an observation `o`.
To maintain consistency, we call everything a state, and assume that an
environment is free to maintain additional (unobserved) internal state.

The `ismdp` query returns true when the environment is MDP, and false otherwise.
"""
ismdp(env::AbstractEnvironment) = false

"""
    maxsteps(env)::Int

Return the max steps in single episode.
Default is `0` (unlimited).
"""
maxsteps(env::AbstractEnvironment) = 0

# ----------------------------------------------------------------
# Implement this interface for a new policy

abstract type AbstractPolicy end

"""
    a = action(policy, r, s, A)

Take in the last reward `r`, current state `s`,
and set of valid actions `A = actions(env, s)`,
then return the next action `a`.

Note that a policy could do a 'sarsa-style' update simply by saving the last state and action `(s,a)`.
"""
function action end

# ----------------------------------------------------------------
# concrete implementations

# include("episodes.jl")
include("episodes/iterators.jl")
include("states.jl")
include("policies/policies.jl")
include("solvers.jl")

include("envs/cartpole.jl")
include("envs/pendulum.jl")
include("envs/mountain_car.jl")
include("envs/multi-armed-bandit.jl")

@reexport using .MultiArmedBanditEnv

# ----------------------------------------------------------------
# a keyboard action space

struct KeyboardAction
  key
end

mutable struct KeyboardActionSet{T} <: AbstractSet{T}
  keys::Vector
end

LearnBase.randtype(s::KeyboardActionSet) = KeyboardAction
Base.rand(s::KeyboardActionSet) = KeyboardAction(rand(s.keys))
Base.in(a::KeyboardAction, s::KeyboardActionSet) = a.key in s.keys
Base.length(s::KeyboardActionSet) = 1

# ----------------------------------------------------------------
# a mouse/pointer action space

struct MouseAction
  x::Int
  y::Int
  button::Int
end

mutable struct MouseActionSet{T} <: AbstractSet{T}
  screen_width::Int
  screen_height::Int
  button::DiscreteSet{Vector{Int}}
end

LearnBase.randtype(s::MouseActionSet) = MouseAction
Base.rand(s::MouseActionSet) =
  MouseAction(rand(1:s.screen_width), rand(1:s.screen_height), rand(s.button))
Base.in(a::MouseAction, s::MouseActionSet) =
  a.x in 1:s.screen_width && a.y in 1:s.screen_height && a.button in s.button
Base.length(s::MouseActionSet) = 1

# ----------------------------------------------------------------

end  # module Reinforce


text:
struct CorrelationMatrixCholesyFactor{K}  <: AbstractParameter end

memory_length(::CorrelationMatrixCholesyFactor{K}) where {K} = (K * (K-1)) >>> 1
memory_length(::Type{<:CorrelationMatrixCholesyFactor{K}}) where {K} = (K * (K-1)) >>> 1


text:
using Revise
using COAST

n_iter = 100

# general params

U238 = 50.0*1e-6
Th232 = 0.0*1e-6

# damage params
c0=6.24534
c1=-0.11977
c2=-314.937
c3=-14.2868
alpha=-0.05721

# functional form, amorphous fraction, mean distance to track parameters
Ea=165.0 # kJ/mol
D0=193188 # cm2/s
D0N17=0.0034 # cm2/s
EaN17=71.0 # kJ/mol
R=COAST.Rjoules/1000 # kJ/mol
Ba=5.48e-19 #g amorphized per alpha event
interconnection=3
SV=1.669 #nm^-1 track surface to volume ratio
lint_lattice=45920.0 #nm, extrapolated to a zircon with 1e14 alphas/g

# preprocess concentrations
(U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol) = conc_to_atoms_per_volume(U238,Th232,density=3.2,U38_35_ratio=137.818)

L = 60e-6
n_T_pts = 221
T2 = collect(LinRange(140.0,40.0,floor(Int,n_T_pts)).+273.15)
T = collect(LinRange(40.0,140.0,ceil(Int,n_T_pts)).+273.15)
T = vcat(T2,T)
times = collect(LinRange(160.0,120.1,ceil(Int,n_T_pts)+1).*3.15576e7*1e6)
times2 = collect(LinRange(120.0,0.0,ceil(Int,n_T_pts)).*3.15576e7*1e6)
times = vcat(times,times2)

he_est = zrdaam_forward_diffusion(n_iter,c0,c1,c2,c3,alpha,Ba,interconnection,SV,lint_lattice,D0,Ea,R,D0N17,EaN17,U238_mol,U235_mol,Th232_mol,L,
  times...,T...)

pre_he = (8*(U238_mol*exp(66.0*1e6*sec_in_yrs/τ38)-U238_mol)+ # HeFTy output
  7*(U235_mol*exp(43.4*1e6*sec_in_yrs/τ35)-U235_mol)+ # HeFTy output
  6*(Th232_mol*exp(43.4*1e6*sec_in_yrs/τ32)-Th232_mol)) # HeFTy output

print("Zirc is: $he_est\n")
print("predicted is: $pre_he\n")

# apatite test
alpha = 0.04672
c0 = 0.39528
c1 = 0.01073
c2 = -65.12969
c3 =  -7.91715
rmr0 = 0.79
eta_q = 0.91
L_dist = 8.1*1e-4 # cm (!)

# general params
U238 = 50.0*1e-6
Th232 = 0.0*1e-6
R = 1.9872*1e-3
n_T_pts = 121
n_iter = 100.0

# preprocess concentrations
(U238_V,U235_V,Th232_V,U238_mol,U235_mol,Th232_mol) = conc_to_atoms_per_volume(U238,Th232,density=3.2,U38_35_ratio=137.818)

"""
check if RDAAM reproduces vanilla Durango
from Farley et al, 00
"""
L1 = 90*1e-4
logD0_a2 = log10(10^1.5/(L1^2)) # cm^2/s
Ea = 32.9

T = collect(LinRange(120.0,0.01,ceil(Int,n_T_pts)).+273.15)
times = collect(LinRange(120.0,0.01,ceil(Int,n_T_pts)+1).*3.1558e7*1e6)

"""
calculate erho
rdaam params in diffn
test fig 2, 2009 RDAAM
"""
Rjoules = 8.314
psi = 1e-13
omega = 1e-22
Etrap = 34*1e3 # J/mol
E_L = 122.3*1e3 # J/mol
L = 60*1e-4
L2 = 60*1e-4
log10D0L_a2_rdaam = log10(exp(9.733)*L^2/L2^2)
T = collect(LinRange(120.0,0.0,ceil(Int,n_T_pts)).+273.15)
times = collect(LinRange(120.0,0.0,ceil(Int,n_T_pts)+1).*3.1558e7*1e6)
mass_he_t2 = rdaam_forward_diffusion(alpha,c0,c1,c2,c3,0.79,eta_q,L_dist,psi,omega,Etrap,Rjoules,E_L,log10D0L_a2_rdaam,n_iter,
               U238_mol,U238_V,U235_mol,U235_V,Th232_mol,Th232_V,L,times...,T...)
pre_time2 = 55.2*1e6
pre_he_t2 = (8*(U238_mol*exp(pre_time2*sec_in_yrs/τ38)-U238_mol)+
             7*(U235_mol*exp(pre_time2*sec_in_yrs/τ35)-U235_mol)+
             6*(Th232_mol*exp(pre_time2*sec_in_yrs/τ32)-Th232_mol))
#print("Ap is: $mass_he_t2\n")
#print("$pre_he_t2\n")


text:
using Test
using clima_demo

@testset "My test" begin
	@test 1 + 1 == 2
end


text:
using transD_GP.GP, Random, PyPlot, Statistics, LinearAlgebra, PositiveFactorizations
##
Random.seed!(12)
f = 1.0/0.5
t = 0:(100*f)^-1:2/f
ntrain = 1
ntries = 4
λ = [0.1]
δ = 0.01
δtry = 5δ
p = 2
demean = true
findmin = true
switchflag = findmin ? -1. : 1.
knownvalue = NaN
acqfun = GP.EI()

y = 100000*switchflag*(2t+ sin.(2*pi*f*t) + sinc.(2*f*(t.-middle(t))))
ynoisy = y + δ*randn(length(t))

tperm = randperm(length(t))
trainidx = tperm[1:ntrain]
ytrain = ynoisy[trainidx]
ttrain = t[trainidx]

ytest, σ2, σ_prior  = GP.GPfit(GP.SqEuclidean(), ytrain', ttrain', t', λ.^2, δtry, p=p, demean=demean)
AF = GP.getAF(acqfun, ytrain, vec(ytest), diag(σ2), findmin=findmin, knownvalue=knownvalue)
f, ax = plt.subplots(2, ntries+1, figsize=(12,5), sharex="col", sharey="row")
ax[1].plot(t,y, "-k")
ax[1].plot(t,ytest)
ax[1].plot(ttrain,ytrain, ".m",markersize=20)
ax[1].plot(t,ytest+2sqrt.(diag(σ2)),"--r",alpha=0.2)
ax[1].plot(t,ytest-2sqrt.(diag(σ2)),"--b",alpha=0.2)
ax[1].grid()
ax[2].plot(t, AF)
ax[2].grid()
ax[1].set_title("step 0")
##
for i = 1:ntries
    global AF
    nextpos = argmax(AF)
    push!(ytrain, ynoisy[nextpos])
    push!(ttrain, t[nextpos])
    ytest, σ2, σ_prior  = GP.GPfit(GP.SqEuclidean(), ytrain', ttrain', t', λ.^2, δtry, p=p, demean=demean)
    AF = GP.getAF(acqfun, ytrain, vec(ytest), diag(σ2), findmin=findmin, knownvalue=knownvalue)
    ax[2i+1].plot(t,y, "-k")
    ax[2i+1].plot(ttrain,ytrain, ".m",markersize=20)
    ax[2i+1].plot(t,ytest)
    ax[2i+1].plot(t,ytest+2sqrt.(diag(σ2)),"--r",alpha=0.2)
    ax[2i+1].plot(t,ytest-2sqrt.(diag(σ2)),"--b",alpha=0.2)
    ax[2i+1].grid()
    ax[2(i+1)].plot(t, AF)
    ax[2(i+1)].grid()
    ax[2i+1].set_title("step $i")
end
# savefig("fig.png",dpi=300)

##
figure()
s1=subplot(211, aspect="auto")
imshow(σ_prior, extent=[t[1],t[end],t[end],t[1]])
s2=subplot(212, sharex=s1, sharey=s1)
imshow(σ2, extent=[t[1],t[end],t[end],t[1]])

figure()
plot(t,y, "-k",label="true")
plot(ttrain,ytrain, "+m", label="train",markersize=10)
nsamp = length(t)
nreal=50
realizations = randn(nreal,nsamp)*cholesky(Positive, var(ytest)*σ2 + 1e-12*Matrix(I,nsamp,nsamp)).U
realizations = ytest .+ realizations'
plot(t,realizations,alpha=0.2)


text:
type BitsContainer <: RoaringContainer
    card::UInt
    data::Vector{UInt64}
end

function BitsContainer()
    Vector{UInt64}()
    return BitsContainer()
end

bitset_container_t *bitset_container_create(void) {
    bitset_container_t *bitset =
        (bitset_container_t *)calloc(1, sizeof(bitset_container_t));

    if (!bitset) {
        return NULL;
    }
    // sizeof(__m256i) == 32
    bitset->array = (uint64_t *) aligned_malloc(32, sizeof(uint64_t) * BITSET_CONTAINER_SIZE_IN_WORDS);
    if (! bitset->array) {
        free(bitset);
        return NULL;
    }
    bitset_container_clear(bitset);
    return bitset;
}


text:
@testset "irreducibles $p" for p = (2, 3)
    Z = ZZ/p
    n = 5
    irrn = collect(irreducibles(Z[:x], n))
    @test irreducible(Z[:x], n) == irrn[1]
    @test length(irrn) == necklace(p, n)
    m = 7
    irrm = collect(irreducibles(Z[:x], m))
    f = prod(irrn[2:4]) * prod(irrm[3:4]) * irrn[2]
    @test isreducible(f)
    ff = factor(f)
    @test length(ff) == 5
    @test prod(ff) == f

    @test length(irrn) == num_irreducibles(Z[:x], n)
    @test length(irrn) == necklace(p, n)
    @test length(irrm) == necklace(p, m)

    redn = collect(reducibles(Z[:x], n))
    @test length(redn) + length(irrn) == p^n
end

@testset "factor non-monic" begin
    Z = ZZ/3
    x = monom(Z[:x])
    p = (-x+1)^9
    @test prod(factor(p)) == p
    p = zero(p)
    @test factor(p) == [p => 1]
    p = one(p)
    @test factor(p) == [p => 1]
    p = x + 1
    @test factor(p) == [p => 1]
    p = -x + 1
    @test factor(p) == [-one(x) => 1, x - 1 => 1]
    p = -(x+1)^3
    @test factor(p) == [-one(x) => 1, x + 1 => 3]
end


text:
using Manopt, Manifolds, Test

@testset "Gradient Descent" begin
    @testset "allocating Circle" begin
        # Test the gradient descent with
        # the distance function squared
        # on S1, such that we can easily also verify exp and log
        M = Circle()
        r = [-π / 2, π / 4, 0.0, π / 4]
        s = r
        F(M, x) = 1 / 10 * sum(distance.(Ref(M), s, Ref(x)) .^ 2)
        gradF(M, x) = 1 / 5 * sum(-log.(Ref(M), Ref(x), s))
        o = gradient_descent!(
            M,
            F,
            gradF,
            s[1];
            stopping_criterion=StopWhenAny(
                StopAfterIteration(200), StopWhenChangeLess(10^-16)
            ),
            stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.99, 0.1),
            debug=[:Iteration, " ", :Cost, :Stop, 100, "\n"],
            record=[:Iteration, :Cost, 1],
            return_options=true,
        )
        x = get_solver_result(o)
        x2 = gradient_descent!(
            M,
            F,
            gradF,
            s[1];
            stopping_criterion=StopWhenAny(
                StopAfterIteration(200), StopWhenChangeLess(10^-16)
            ),
            stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.99, 0.1),
        )
        @test x == x2
        x3 = gradient_descent!(
            M,
            F,
            gradF,
            s[1];
            stopping_criterion=StopWhenAny(
                StopAfterIteration(1000), StopWhenChangeLess(10^-16)
            ),
            stepsize=NonmonotoneLinesearch(
                1.0,
                ExponentialRetraction(),
                ParallelTransport(),
                0.99,
                0.1,
                2,
                1e-7,
                π / 2,
                :direct,
            ),
            debug=[:Stop],
        )
        @test isapprox(x, x3; atol=1e-13)
        x4 = gradient_descent!(
            M,
            F,
            gradF,
            s[1];
            stopping_criterion=StopWhenAny(
                StopAfterIteration(1000), StopWhenChangeLess(10^-16)
            ),
            stepsize=NonmonotoneLinesearch(
                1.0,
                ExponentialRetraction(),
                ParallelTransport(),
                0.99,
                0.1,
                2,
                1e-7,
                π / 2,
                :inverse,
            ),
            debug=[:Stop],
        )
        @test isapprox(x, x4; atol=1e-13)
        x5 = gradient_descent!(
            M,
            F,
            gradF,
            s[1];
            stopping_criterion=StopWhenAny(
                StopAfterIteration(1000), StopWhenChangeLess(10^-16)
            ),
            stepsize=NonmonotoneLinesearch(
                1.0,
                ExponentialRetraction(),
                ParallelTransport(),
                0.99,
                0.1,
                2,
                1e-7,
                π / 2,
                :alternating,
            ),
            debug=[:Stop],
        )
        @test isapprox(x, x5; atol=1e-13)
        x6 = gradient_descent!(
            M,
            F,
            gradF,
            s[1];
            stopping_criterion=StopWhenAny(
                StopAfterIteration(1000), StopWhenChangeLess(10^-16)
            ),
            direction=Nesterov(s[1]),
            debug=[:Stop],
        )
        @test isapprox(x, x6; atol=1e-13)

        @test_logs (
            :warn,
            string(
                "The strategy 'indirect' is not defined. The 'direct' strategy is used instead.",
            ),
        ) NonmonotoneLinesearch(M; strategy=:indirect)
        @test_throws DomainError NonmonotoneLinesearch(M; min_stepsize=0.0)
        @test_throws DomainError NonmonotoneLinesearch(
            M; min_stepsize=π / 2, max_stepsize=π / 4
        )
        @test_throws DomainError NonmonotoneLinesearch(M; memory_size=0)

        rec = get_record(o)
        # after one step for local enough data -> equal to real valued data
        @test abs(x - sum(r) / length(r)) ≈ 0 atol = 5 * 10.0^(-14)
        # Test Fallbacks -> we can't do steps with the wrong combination
        p = SubGradientProblem(M, F, gradF)
        o = GradientDescentOptions(
            s[1]; stopping_criterion=StopAfterIteration(20), stepsize=ConstantStepsize(M)
        )
        @test_throws MethodError initialize_solver!(p, o)
        @test_throws MethodError step_solver!(p, o, 1)
    end
    @testset "mutating Sphere" begin
        M = Sphere(2)
        north = [0.0, 0.0, 1.0]
        pre_pts = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], [0.0, -1.0, 0.0]]
        pts = exp.(Ref(M), Ref(north), pre_pts)
        F(M, x) = 1 / 8 * sum(distance.(Ref(M), pts, Ref(x)) .^ 2)
        gradF(M, x) = 1 / 4 * sum(-log.(Ref(M), Ref(x), pts))
        n2 = gradient_descent(M, F, gradF, pts[1])
        @test !isapprox(M, pts[1], n2) # n2 is newly allocated and not pts[1]
        @test isapprox(M, north, n2)
        n3 = gradient_descent(
            M, F, gradF, pts[1]; direction=MomentumGradient(M, pts[1]), debug=[]
        )
        @test isapprox(M, north, n3)
        n4 = gradient_descent(M, F, gradF, pts[1]; direction=AverageGradient(M, pts[1], 5))
        @test isapprox(M, north, n4; atol=1e-7)
    end
    @testset "Warning when cost increases" begin
        M = Sphere(2)
        q = 1 / sqrt(2) .* [1.0, 1.0, 0.0]
        F(M, p) = distance(M, p, q) .^ 2
        # chosse a wrong gradient such that ConstantStepsize yields an increase
        gradF(M, p) = -grad_distance(M, q, p)
        # issues three warnings
        @test_logs (:warn,) (:warn,) (:warn,) gradient_descent(
            M,
            F,
            gradF,
            1 / sqrt(2) .* [1.0, -1.0, 0.0];
            debug=[DebugWarnIfCostIncreases(:Once)],
        )
    end
end


text:
using HorizonSideRobots

function obhod(r::Robot,side::HorizonSide)
    k=0
    j=0
    while isborder(r,side)==true && isborder(r,next(side))==false
        move!(r,next(side))
        k+=1
    end
    if isborder(r,side)==true && isborder(r,next(side))==true
        for _ in 1:k
            move!(r,inverse(next(side)))
        end
        k=0
    end
    while isborder(r,side)==true && isborder(r,inverse(next(side)))==false
        move!(r,inverse(next(side)))
        j+=1
    end
    if isborder(r,side)==true && isborder(r,inverse(next(side)))==true
        for _ in 1:j
            move!(r,next(side))
        end
    end
    if k==0
        k=j
        side1 = next(side)
    else
        side1 = inverse(next(side))
    end
    if isborder(r,side)==false
        print(side)
        move!(r,side)
        while (isborder(r,next(inverse(side)))==true &&(isborder(r,next(side))==false || isborder(r,inverse(next(side)))==false))
            move!(r,side)
        end
        for _ in 1:k
            move!(r,side1)
        end
    end
end

function find_mrk_spiral(r)

counter = 1
while ismarker(r)==false

    for i in 1:counter
        obhod(r, Ost)
        if ismarker(r)==true
            break
        end
    end
    if ismarker(r)==true
        break
    end

    for i in 1:counter
        obhod(r, Nord)
        if ismarker(r)==true
            break
        end
    end
    if ismarker(r)==true
        break
    end

    counter += 1

    for i in 1:counter
        obhod(r, West)
        if ismarker(r)==true
            break
        end
    end
    if ismarker(r)==true
        break
    end

    for i in 1:counter
        obhod(r, Sud)
        if ismarker(r)==true
            break
        end
    end
    if ismarker(r)==true
        break
    end

    counter += 1

end
end

find_mrk_spiral(r)


text:
"""
    MPO{L,T<:Number}

Matrix product operator for `L` sites with numerical entries of type `T`.

The tensors are stored with the indices `[a, m, n, b]` corresponding to

       n
       |
    a -#- b
       |
       m

where the physical indices (`m`, `n`) are between the bond indices (`a`, `b`).
The end tensors also have this structure, so they must be capped explicitly.
"""
struct MPO{L,T<:Number}
    "Tensors making up the MPO."
    tnsrs::Vector{Array{T,4}}
end

"""
    MPO(x::Array{T,4}, L)

Create an `MPO` for `L` sites with all interior sites containing the tensor
`x`. The tensor is assumed to have the usual matrix-of-operators structure,
with the first two indices being the bond (matrix) dimension and the last two
indices being the physical (operator) dimension. The first and last sites only
use the second (not last!) row and first column of `x`, respectively.
"""
function MPO(x::Array{T,4}, L) where {T}
    L >= 2 || throw(DomainError(L, "At least 2 sites."))

    tnsrs = Vector{Array{T,4}}(undef, L)
    # Row vector.
    tnsrs[1] = permutedims(x[2:2, :, :, :], (1, 3, 4, 2))
    for i in 2:(L-1)
        # Matrix.
        tnsrs[i] = permutedims(x, (1, 3, 4, 2))
    end
    # Column vector.
    tnsrs[L] = permutedims(x[:, 1:1, :, :], (1, 3, 4, 2))

    MPO{L,T}(tnsrs)
end

"""
    compress!(o::MPO, cutoff_max::Float64)

Compress `o` by throwing away singular values below `cutoff_max` at each bond.
"""
function compress!(o::MPO{L}, cutoff_max::Float64) where {L}
    bond_dimensions = Array{Int}(undef, L-1)

    for range in [1:(L-2), (L-1):-1:1]
        for site in range
            # Combine adjacent tensors:
            #
            #        n     p
            #        |     |
            #     a -#- b -#- c
            #        |     |
            #        m     o
            @tensor M[a, m, n, o, p, c] :=
                o.tnsrs[site][a, m, n, b] * o.tnsrs[site+1][b, o, p, c]
            M_mat = reshape(M, size(M, 1) * size(M, 2) * size(M, 3),
                               size(M, 4) * size(M, 5) * size(M, 6))

            U, S, V = svd(M_mat)
            # Operators are not normalized in general, so we should not expect
            # that sum(S.^2) == 1.
            trunc_total = sum(S.^2)
            trunc_len = length(S)
            trunc_cutoff = 0.0
            while (1 < trunc_len &&
                   trunc_cutoff + S[trunc_len]^2/trunc_total <= cutoff_max)
                trunc_cutoff += S[trunc_len]^2/trunc_total
                trunc_len -= 1
            end
            A = U[:, 1:trunc_len]
            B = diagm(S[1:trunc_len]) * V[:, 1:trunc_len]'
            o.tnsrs[site] = reshape(A, size(M, 1), size(M, 2), size(M, 3),
                                       trunc_len)
            o.tnsrs[site+1] = reshape(B, trunc_len,
                                         size(M, 4), size(M, 5), size(M, 6))

            bond_dimensions[site] = trunc_len
        end
    end

    mean(bond_dimensions)
end


text:
using .Nodes
using ..Partitioners

# TODO: Allow the instantiation without the specification of width (Could be
#       determined from the training data)
struct Discriminator{P <: AbstractPartitioner,T <: AbstractNode} <: AbstractModel
    partitioner::P
    nodes::Vector{T}
end

function Discriminator{P,T}(partitioner::P; kargs...) where {P <: AbstractPartitioner,T <: AbstractNode}
    nodes = [T(;kargs...) for _ in 1:length(partitioner)]

    Discriminator{P,T}(partitioner, nodes)
end

function Discriminator{RandomPartitioner,T}(width::Int, n::Int; seed::Union{Nothing,Int}=nothing, kargs...) where {T <: AbstractNode}
    partitioner = RandomPartitioner(width, n; seed)

    Discriminator{RandomPartitioner,T}(partitioner; kargs...)
end

function Discriminator{P,T}(X::U, partitioner::P; kargs...) where {P <: AbstractPartitioner,T <: AbstractNode,U <: AbstractVecOrMat{Bool}}
    d = Discriminator{P,T}(partitioner; kargs...)

    train!(d, X)

    return d
end

function Discriminator{RandomPartitioner,T}(X::U, n::Int; seed::Union{Nothing,Int}=nothing, kargs...) where {T <: AbstractNode,U <: AbstractVecOrMat{Bool}}
    d = Discriminator{RandomPartitioner,T}(size(X)[end], n; seed, kargs...)

    train!(d, X)

    return d
end

## Aliases and convenience constructors
# Default node is DictNode with random partitioning
Discriminator(args...; kargs...) = Discriminator{RandomPartitioner,DictNode}(args...; kargs...)

const BitDiscriminator        = Discriminator{RandomPartitioner,DictNode{BitVector}}
const BleachingDiscriminator  = Discriminator{RandomPartitioner,AccNode} # This should be the default classification discriminator

RegressionDiscriminator = Discriminator{RandomPartitioner,RegressionNode}
GeneralizedRegressionDiscriminator = Discriminator{RandomPartitioner,GeneralizedRegressionNode}

# Generic functions
function train!(d::Discriminator{<:AbstractPartitioner,<:AbstractClassificationNode}, X::T) where {T <: AbstractVecOrMat{Bool}}
    for (node, x) in Iterators.zip(d.nodes, partition(d.partitioner, X))
        Nodes.train!(node, x)
    end

    return nothing
end

function train!(d::Discriminator{<:AbstractPartitioner,<:AbstractRegressionNode}, X, y)
    for (node, x) in Iterators.zip(d.nodes, partition(d.partitioner, X))
        Nodes.train!(node, x, y)
    end

    return nothing
end

initial_response(::AbstractVector{Bool}) = zero(Int)
initial_response(X::AbstractMatrix{Bool}) = zeros(Int, size(X, 1))

# TODO: Output response in relative terms? (i.e. divide by the number of nodes)
function predict(d::Discriminator{<:AbstractPartitioner,<:AbstractClassificationNode}, X::AbstractVecOrMat{Bool}; b::Int=0)
    response = initial_response(X)

    for (node, x) in Iterators.zip(d.nodes, partition(d.partitioner, X))
        response += predict(node, x; b)
    end

    return response
end

function predict(d::Discriminator{<:AbstractPartitioner,<:AbstractRegressionNode}, X::AbstractMatrix{Bool})
    [predict(d, x) for x in eachrow(X)]
end

# TODO: In the context of discriminators, predict and predict_response are the same thing.
#       Implement predict_response as an alias to predict

################################################################################
# Specialized training and prediction methods for the regresion variant

function predict(d::Discriminator{P,RegressionNode}, X::AbstractVector{Bool}) where {P <: AbstractPartitioner}
    partial_count = zero(Int)
    estimate = zero(Float64)

    for (node, x) in Iterators.zip(d.nodes, partition(d.partitioner, X))
        count, sum = predict(node, x)

        if node.γ != 1.0
            count = (1 - node.γ^count) / (1 - node.γ)
        end

        partial_count += count
        estimate += sum
    end

    return partial_count == 0 ? estimate : estimate / partial_count
end

################################################################################
# Specialized training and prediction methods for the generalized regresion discriminator

# This expects that α is a constant
function predict(d::Discriminator{P,GeneralizedRegressionNode}, X::AbstractVector{Bool}) where {P <: AbstractPartitioner}
    running_numerator = zero(Float64)
    running_denominator = zero(Float64)

    for (node, x) in Iterators.zip(d.nodes, partition(d.partitioner, X))
        count, estimate = predict(node, x)

        α = node.α
        if count != 0
            running_numerator += estimate / α
            running_denominator += 1 / α
        end
    end

    return running_denominator == 0 ? zero(Float64) : running_numerator / running_denominator
end

max_stars_repo_name: RafaelFK/ramnet

max_issues_repo_name: RafaelFK/ramnet

max_forks_repo_name: RafaelFK/ramnet