text:
# auxiliary functions for all the other methods
# originally in functions.jl
# Claudia February 2015
#####################

function setCHECKNET(b::Bool)
    global CHECKNET
    CHECKNET = b
    CHECKNET && warn("PhyloNetworks.CHECKNET is true: will slow snaq! down.")
    b || println("PhyloNetworks.CHECKNET set to false")
end

# ----- aux general functions ---------------

#based in coupon's collector: E+sqrt(V)
function coupon(n::Number)
    return n*log(n) + n
end

function binom(n::Number,k::Number)
    n >= k || return 0
    n == 1 && return 1
    k == 0 && return 1
    binom(n-1,k-1) + binom(n-1,k) #recursive call
end

function approxEq(a::Number,b::Number,absTol::Number,relTol::Number)
    if(a<eps() || b<eps())
        abs(a-b) < absTol
    else
        abs(a-b) < relTol*eps(abs(a)+abs(b))
    end
end

approxEq(a::Number,b::Number) = approxEq(a,b,1e-5,100)

# isEqual functions: to test if 2 edges (or 2 nodes etc.) "look" alike.
#                    Useful after a deepcopy of a network.
# For nodes (or edges etc.) in the same network, use instead n1 == n2 or n1 != n2.
function isEqual(n1::Node,n2::Node)
    return (n1.number == n2.number && approxEq(n1.gammaz,n2.gammaz) && n1.inCycle == n2.inCycle)
end

function isEqual(n1::Edge,n2::Edge)
    return (n1.number == n2.number && approxEq(n1.length,n2.length))
end

function isEqual(net1::HybridNetwork, net2::HybridNetwork)
    result = true
    result &= (net1.numTaxa == net2.numTaxa)
    result &= (net1.numNodes == net2.numNodes)
    result &= (net1.numEdges == net2.numEdges)
    ## result &= (net1.node == net2.node)
    ## result &= (net1.edge == net2.edge)
    result &= (net1.root == net2.root)
    result &= (net1.names == net2.names)
##    result &= (net1.hybrid == net2.hybrid)
    result &= (net1.numHybrids == net2.numHybrids)
##    result &= (net1.leaf == net2.leaf)
    result &= (net1.ht == net2.ht)
    result &= (net1.numht == net2.numht)
    result &= (net1.numBad == net2.numBad)
    result &= (net1.hasVeryBadTriangle == net2.hasVeryBadTriangle)
    result &= (net1.index == net2.index)
    result &= (net1.loglik == net2.loglik)
    return result
end

#------------- functions to allow for ------------#
#              missing values (lengths or gammas) #

# adds x+y but interprets -1.0 as missing: so -1.0 + x = -1.0 here.
function addBL(x::Number,y::Number)
    (x==-1.0 || y==-1.0) ? -1.0 : x+y
end

#------------- EDGE functions --------------------#

# warning: node needs to be defined as hybrid before adding to a
#          hybrid edge. First, an edge is defined as hybrid, and then
#          the nodes are added to it. If the node added is leaf, the
#          edge length is set unidentifiable (as it is external edge)
function setNode!(edge::Edge, node::Node)
    global DEBUG
    size(edge.node,1)  !=  2 || error("vector of nodes already has 2 values");
    push!(edge.node,node);
    if(size(edge.node,1) == 1)
        if(edge.hybrid)
            if(node.hybrid)
                edge.isChild1 = true;
            else
                edge.isChild1 = false;
	    end
        end
        if(node.leaf)
            edge.istIdentifiable = false;
        else
            edge.istIdentifiable = true
        end
    else
        if(node.leaf)
            !edge.node[1].leaf || error("edge $(edge.number) has two leaves")
            edge.istIdentifiable = false;
        else
            if(edge.hybrid)
	        if(node.hybrid)
                    if(DEBUG)
                        !edge.node[1].hybrid || println("hybrid edge $(edge.number) has two hybrid nodes");
                    end
                    edge.isChild1 = false;
	        else
	            edge.node[1].hybrid || error("hybrid edge $(edge.number) has no hybrid nodes");
	            edge.isChild1 = true;
	        end
            else #edge is tree
                if(!edge.node[1].leaf)
                    if(!node.hybrid && !edge.node[1].hybrid)
                        if(edge.fromBadDiamondI)
                            edge.istIdentifiable = false
                        else
                            edge.istIdentifiable = true
                        end
                    else
                        if(node.hybrid && (node.isBadDiamondI || node.isBadDiamondII || node.isBadTriangle))
                            edge.istIdentifiable = false
                        elseif(edge.node[1].hybrid && (edge.node[1].isBadDiamondI ||edge.node[1].isBadDiamondII || edge.node[1].isBadTriangle))
                            edge.istIdentifiable = false
                        else
                            edge.istIdentifiable = true
                        end
                    end
                else
                    edge.istIdentifiable = false
                end
            end
        end
    end
end

# warning: node needs to be defined as hybrid before adding to a hybrid edge.
#          First, an edge is defined as hybrid, and then the nodes are added to it.
#          If there is a leaf in node, the edge.istIdentifiable=false
function setNode!(edge::Edge,node::Array{Node,1})
    size(node,1) ==  2 || error("vector of nodes must have exactly 2 values")
    edge.node = node;
    if(edge.hybrid)
      if(node[1].hybrid)
          edge.isChild1 = true;
      else
          node[2].hybrid || error("hybrid edge without hybrid node");
          edge.isChild1 = false;
      end
    end
    if(edge.node[1].leaf || edge.node[2].leaf)
        edge.istIdentifiable = false;
    else
        edge.istIdentifiable = true;
    end
end

"""
    getChild(edge::Edge)

Return child node using the `isChild1` attribute of the edge.
"""
@inline function getChild(edge::Edge)
    edge.node[edge.isChild1 ? 1 : 2]
end

"""
    getParent(edge::Edge)

Return parent node using the `isChild1` attribute of the edge.
"""
@inline function getParent(edge::Edge)
    edge.node[edge.isChild1 ? 2 : 1]
end

"""
    getPartner(edge::Edge)

Return hybrid partner of edge, that is, hybrid edge pointing
to the same child as `edge`. Assumes correct `isChild1` attributes.
Assumes no in-coming polytomy: a node has 0, 1 or 2 parents, no more.
"""
@inline function getPartner(edge)
    node = getChild(edge)
    getPartner(edge, node)
end
@inline function getPartner(edge::Edge, node::Node)
    for e in node.edge
        if e.hybrid && e != edge && node == getChild(e)
            return e
        end
    end
    error("did not find a partner for edge $(edge.number)")
end
# -------------- NODE -------------------------#

function setEdge!(node::Node,edge::Edge)
   push!(node.edge,edge);
   all((e->!e.hybrid), node.edge) ? node.hasHybEdge = false : node.hasHybEdge = true;
end

function getOtherNode(edge::Edge,node::Node)
  isequal(edge.node[1],node) ? edge.node[2] : edge.node[1]
end

"""
    getMajorParent(node::Node)
    getMinorParent(node::Node)

Return major or minor parent of a node using the `isChild1` field of edges
(and assuming correct `isMajor` field).
"""
@inline function getMajorParent(node::Node)
    for e in node.edge
        if e.isMajor && node == getChild(e)
            return getParent(e)
        end
    end
    error("could not find major parent of node $(node.number)")
end
@doc (@doc getMajorParent) getMinorParent
@inline function getMinorParent(node::Node)
    for e in node.edge
        if !e.isMajor && node == getChild(e)
            return getParent(e)
        end
    end
    error("could not find minor parent of node $(node.number)")
end

# -------------- NETWORK ----------------------- #

function getIndex(node::Node, net::Network)
    i = 1;
    while(i<= size(net.node,1) && !isEqual(node,net.node[i]))
        i = i+1;
    end
    i>size(net.node,1)?error("node $(node.number) not in network"):return i;
end

function getIndex(edge::Edge, net::Network)
    i = 1;
    while(i<= size(net.edge,1) && !isEqual(edge,net.edge[i]))
        i = i+1;
    end
    i>size(net.edge,1)?error("edge $(edge.number) not in network"):return i;
end

function getIndex(edge::Edge, edges::Vector{Edge})
    i = 1;
    while(i<= size(edges,1) && !isEqual(edge,edges[i]))
        i = i+1;
    end
    i>size(edges,1)?error("edge $(edge.number) not in array of edges"):return i;
end

function getIndex(bool::Bool, array::Array{Bool,1})
    i = 1;
    while(i<= size(array,1) && !isequal(bool,array[i]))
        i = i+1;
    end
    i>size(array,1)?error("$(bool) not in array"):return i;
end

## function getIndex(bool::Bool, array::Array{Bool,1})
##     i = 1;
##     while(i<= size(array,1) && !isequal(bool,array[i]))
##         i = i+1;
##     end
##     i>size(array,1)?error("$(bool) not in array"):return i;
## end

function getIndex(bool::Bool, array::Array{Any,1})
    i = 1;
    while(i<= size(array,1) && !isequal(bool,array[i]))
        i = i+1;
    end
    i>size(array,1)?error("$(bool) not in array"):return i;
end

# aux function to find the index of a string in a
# string array
function getIndex(name::AbstractString, array::Array{String,1})
    i = 1;
    while(i<= size(array,1) && !isequal(name,array[i]))
        i = i+1;
    end
    i>size(array,1)?error("$(name) not in array"):return i;
end

# aux function to find the index of a int in an int array.
# But findfirst can do that as well, and probably more efficiently (returning 0 if not found)
function getIndex(name::Integer, array::Array{Int,1})
    i = 1;
    while(i<= size(array,1) && !isequal(name,array[i]))
        i = i+1;
    end
    i>size(array,1)?error("$(name) not in array"):return i;
end

# aux function to find the index of a node in a
# node array
function getIndex(name::Node, array::Array{Node,1})
    i = 1;
    while(i<= size(array,1) && !isequal(name,array[i]))
        i = i+1;
    end
    i>size(array,1)?error("$(name.number) not in array"):return i;
end

function getIndexNode(number::Integer,net::Network)
    ind = findfirst([number==n.number for n in net.node])
    if ind==0
        error("node number not in net.node")
    end
    return ind
end

function getIndexEdge(number::Integer,net::Network)
    ind = findfirst([number==n.number for n in net.edge])
    if ind==0
        error("edge number not in net.edge")
    end
    return ind
end

# find the index of an edge in node.edge
function getIndexEdge(edge::Edge,node::Node)
    getIndex(true,[isequal(edge,e) for e in node.edge])
end

# find the index of an edge with given number in node.edge
function getIndexEdge(number::Integer,node::Node)
    getIndex(true,[isequal(edge,e) for e in node.edge])
end

# find the index of a node in edge.node
function getIndexNode(edge::Edge,node::Node)
    size(edge.node,1) == 2 || warn("this edge $(edge.number) has more or less than 2 nodes: $([n.number for n in edge.node])")
    if(isequal(node,edge.node[1]))
        return 1
    elseif(isequal(node,edge.node[2]))
        return 2
    else
        error("node not in edge.node")
    end
end

# function to find hybrid index in net.hybrid
function getIndexHybrid(node::Node, net::Network)
    node.hybrid || error("node $(node.number) is not hybrid so it cannot be in net.hybrid")
    i = 1;
    while(i<= size(net.hybrid,1) && !isEqual(node,net.hybrid[i]))
        i = i+1;
    end
    i>size(net.hybrid,1)?error("hybrid node not in network"):return i;
end

# function to find leaf index in qnet.leaf
function getIndexLeaf(node::Node, net::Network)
    if(node.leaf)
        i = 1;
        while(i<= size(net.leaf,1) && !isEqual(node,net.leaf[i]))
            i = i+1;
        end
        i>size(net.leaf,1)?error("leaf node not in network"):return i;
    else
        error("node $(node.number) is not leaf so it cannot be in net.leaf")
    end
end

# function that given a hybrid node, it gives you the minor hybrid edge
function getHybridEdge(node::Node)
    node.hybrid || error("node $(node.number) is not hybrid node, cannot get hybrid edges")
    a = nothing;
    for e in node.edge
        (e.hybrid && !e.isMajor) ? a = e : nothing;
    end
    isa(a,Void) ? error("hybrid node $(node.number) does not have minor hybrid edge, edges: $([e.number for e in node.edge])") : return a
end

# function that given two nodes, it gives you the edge that connects them
# returns error if they are not connected by an edge
function getConnectingEdge(node1::Node,node2::Node)
    found = false;
    i = 1;
    while(i<= size(node1.edge,1) && !found)
        if(isequal(getOtherNode(node1.edge[i],node1),node2))
            found = true;
        end
        i = i+1;
    end
    if(found)
        return node1.edge[i-1]
    else
        error("nodes not connected")
    end
end

# function to check in an edge is in an array by comparing
# the edges numbers (uses isEqual)
# needed for updateHasEdge
function isEdgeNumIn(edge::Edge,array::Array{Edge,1})
    return all((e->!isEqual(edge,e)), array) ? false : true
end

# function to check in a leaf is in an array by comparing
# the numbers (uses isEqual)
# needed for updateHasEdge
function isNodeNumIn(node::Node,array::Array{Node,1})
    return all((e->!isEqual(node,e)), array) ? false : true
end

# function to push a Node in net.node and
# update numNodes and numTaxa
function pushNode!(net::Network, n::Node)
    push!(net.node,n);
    net.numNodes += 1;
    if(n.leaf)
        net.numTaxa += 1
        push!(net.leaf,n);
    end
    if(n.hybrid)
        pushHybrid!(net,n)
    end
end

# function to push an Edge in net.edge and
# update numEdges
function pushEdge!(net::Network, e::Edge)
    push!(net.edge,e);
    net.numEdges += 1;
end

# function to push a hybrid Node in net.hybrid and
# update numHybrids
function pushHybrid!(net::Network, n::Node)
    if(n.hybrid)
        push!(net.hybrid,n);
        net.numHybrids += 1;
    else
        error("node $(n.number) is not hybrid, so cannot be pushed in net.hybrid")
    end
end

"""
`deleteNode!(net::HybridNetwork, n::Node)`

deletes a Node from a network, i.e. removes it from
net.node, and from net.hybrid or net.leaf as appropriate.
Updates attributes numNodes, numTaxa, numHybrids
(it does not update net.names though).

Warning: if the root is deleted, the new root is arbitrarily set to the
first node in the list. This is intentional to save time because this function
is used frequently in snaq!, which handles semi-directed (unrooted) networks.
"""
function deleteNode!(net::HybridNetwork, n::Node)
    index = 0
    try
        index = getIndex(n,net);
    catch
        error("Node $(n.number) not in network");
    end
    # println("deleting node $(n.number) from net, index $(index).")
    deleteat!(net.node,index);
    net.numNodes -= 1;
    if(net.root == index) # do not check containRoot to save time in snaq!
        net.root = 1      # arbitrary
    elseif(net.root > index)
        net.root -= 1
    end
    if(n.hybrid)
       removeHybrid!(net,n)
    end
    if(n.leaf)
        removeLeaf!(net,n)
    end
end

# function to delete a Node in net.node and
# update numNodes and numTaxa for QuartetNetwork
# if hybrid node, it deletes also from net.hybrid
# and updates numHybrids
# note that net.names is never updated to keep it
# accurate
# if n is leaf, we delete from qnet.leaf
function deleteNode!(net::QuartetNetwork, n::Node)
    index=0
    try
        index = getIndex(n,net);
    catch
        error("Node $(n.number) not in network");
    end
    #println("deleting node $(n.number) from net")
    deleteat!(net.node,index);
    net.numNodes -= 1;
    net.numTaxa -= n.leaf ? 1 : 0;
    if(n.hybrid)
       removeHybrid!(net,n)
    end
    if(n.leaf)
        index = getIndexLeaf(n,net)
        deleteat!(net.leaf,index)
    end
end

# function to delete an Edge in net.edge and
# update numEdges from a HybridNetwork
# added part boolean, default true to check the partition only when part=true
function deleteEdge!(net::HybridNetwork, e::Edge; part=true::Bool)
    if(part)
        if(e.inCycle == -1 && !e.hybrid && !isempty(net.partition) && !isTree(net))
            ind = whichPartition(net,e)
            indE = getIndex(e,net.partition[ind].edges)
            deleteat!(net.partition[ind].edges,indE)
        end
    end
    index = 0
    try
        index = getIndex(e,net);
    catch
        error("Edge $(e.number) not in network");
    end
    #println("delete edge $(e.number) from net")
    deleteat!(net.edge,index);
    net.numEdges -= 1;
end

# function to delete an Edge in net.edge and
# update numEdges from a QuartetNetwork
function deleteEdge!(net::QuartetNetwork, e::Edge)
    index=0
    try
        index = getIndex(e,net);
    catch
        error("Edge not in network");
    end
    #println("delete edge $(e.number) from net")
    deleteat!(net.edge,index);
    net.numEdges -= 1;
end

# function to delete a hybrid Node in net.hybrid and
# update numHybrid
# used when you do not want to delete the actual node
# only remove it from net.hybrid
function removeHybrid!(net::Network, n::Node)
    n.hybrid || error("cannot delete node $(n.number) from net.hybrid because it is not hybrid")
    index = 0
    try
        index = getIndexHybrid(n,net);
    catch
        error("Hybrid Node $(n.number) not in network");
    end
    deleteat!(net.hybrid,index);
    net.numHybrids -= 1;
end

# function to delete a leaf node in net.leaf
# and update numTaxa
function removeLeaf!(net::Network,n::Node)
    if(n.leaf)
        index = 0
        try
            index = getIndexLeaf(n,net)
        catch
            error("Leaf node $(n.number) not in network")
        end
        deleteat!(net.leaf,index)
        net.numTaxa -= 1
    else
        error("cannot delete node $(n.number) from net.leaf because it is not leaf")
    end
end

# function to delete an internal node with only 2 edges
function deleteIntNode!(net::Network, n::Node)
    size(n.edge,1) == 2 || error("node $(n.number) does not have only two edges")
#    isEqual(n,net.node[net.root]) && println("deleting the root $(n.number) because it has only two edges attached")
    index = n.edge[1].number < n.edge[2].number ? 1 : 2;
    edge1 = n.edge[index];
    edge2 = n.edge[index==1?2:1];
    if(!edge1.hybrid && !edge2.hybrid)
        node1 = getOtherNode(edge1,n);
        node2 = getOtherNode(edge2,n);
        removeEdge!(node2,edge2);
        removeNode!(n,edge1);
        setEdge!(node2,edge1);
        setNode!(edge1,node2);
        deleteNode!(net,n);
        deleteEdge!(net,edge2);
    else
        warn("the two edges $([edge1.number,edge2.number]) attached to node $(n.number) must be tree edges to delete node")
        if(edge1.hybrid)
            hybedge = edge1
            otheredge = edge2
        elseif(edge2.hybrid)
            hybedge = edge2
            otheredge = edge1
        end
        othernode = getOtherNode(otheredge,n)
        removeNode!(n,hybedge)
        removeEdge!(othernode,otheredge)
        setEdge!(othernode,hybedge)
        setNode!(hybedge,othernode)
        deleteNode!(net,n)
        deleteEdge!(net,otheredge)
    end
end

# search the hybrid node(s) in network: returns the hybrid node(s)
# in an array
# throws error if no hybrid in network
function searchHybridNode(net::Network)
    a = [n for n in net.node if n.hybrid]
    suma = length(a)
    suma != 0 || error("network has no hybrid node")
    return a
end

# search and returns the hybrid edges in network
# throws error if no hybrid in network
function searchHybridEdge(net::Network)
    a = [n for n in net.edge if n.hybrid]
    suma = length(a)
    suma != 0 || error("network has no hybrid edge")
    return a
end

# print for every edge, nodes, inCycle, containRoot, istIdentifiable
function printEdges(net::QuartetNetwork)
    println("Edge\tNode1\tNode2\tInCycle\tcontainRoot\tistIdentifiable\tLength\tisHybrid\tGamma")
    for e in net.edge
        println("$(e.number)\t$(e.node[1].number)\t$(e.node[2].number)\t$(e.inCycle)\t$(e.containRoot)\t\t$(e.istIdentifiable)\t\t$(round(e.length,2))\t$(e.hybrid)\t$(round(e.gamma,4))")
    end
end

"""
`printEdges(net::HybridNetwork)`

prints the information on the edges of net: edge number, node numbers of nodes attached to it, in which cycle it is contained (-1 if no cycle), can it contain root, is it an identifiable edge, length, is it hybrid, gamma value
"""

function printEdges(net::HybridNetwork)
    if(net.numBad > 0)
        println("net has $(net.numBad) bad diamond I. Some γ and edge lengths t are not identifiable, although their γ * (1-exp(-t)) are.")
    end
    miss = "NA"
    println("Edge\tNode1\tNode2\tInCycle\tcontainRoot\tistIdentitiable\tLength\tisHybrid\tGamma\tisMajor")
    for e in net.edge
        println("$(e.number)\t$(e.node[1].number)\t$(e.node[2].number)\t$(e.inCycle)\t$(e.containRoot)\t\t$(e.istIdentifiable)\t\t$(e.length==-1? miss :round(e.length,2))\t$(e.hybrid)\t\t$(e.gamma == 1.0 && e.hybrid ? miss : round(e.gamma,4))\t$(e.isMajor)")
    end
end

# print for every node, inCycle and edges
"""
`printNodes(net::HybridNetwork)`

prints information on the nodes of net: node number, in which cycle it is contained (-1 if no cycle), is it hybrid, does it has hybrid edges, edges number attached to it
"""
function printNodes(net::Network)
    println("Node\tIn Cycle\tisHybrid\thasHybEdge\tNode label\tisLeaf\tEdges numbers")
    for n in net.node
        print("$(n.number)\t$(n.inCycle)\t\t$(n.hybrid)\t\t$(n.hasHybEdge)\t\t$(n.name)\t\t$(n.leaf)")
        for e in n.edge
            print("\t$(e.number)")
        end
        print("\n")
    end
end

# find the edges for a given hybrid node
# in the order: hybrid major, hybrid minor, tree edge
# if node is tree node with hybrid edges, it returns
# hybrid edge, tree edge in cycle, tree edge not in cycle
# warning: assumes any tree node with hybrid edge has two tree edges
#          one in cycle, the other not in cycle
function hybridEdges(node::Node)
    size(node.edge,1) == 3 || error("node $(node.number) has $(size(node.edge,1)) edges instead of 3");
    if(node.hybrid)
        hybmajor = nothing;
        hybminor = nothing;
        tree = nothing;
        for e in node.edge
            (e.hybrid && e.isMajor) ? hybmajor = e : nothing
            (e.hybrid && !e.isMajor) ? hybminor = e : nothing
            !e.hybrid ? tree = e : nothing
        end
        return hybmajor, hybminor, tree
    elseif(node.hasHybEdge)
        hybrid = nothing;
        treecycle = nothing;
        tree = nothing;
        for e in node.edge
            (e.hybrid) ? hybrid = e : nothing
            (!e.hybrid && e.inCycle != -1) ? treecycle = e : nothing
            (!e.hybrid && e.inCycle == -1) ? tree = e : nothing
        end
        return hybrid, treecycle, tree
    else
        #warn("node $(node.number) is not hybrid $(node.hybrid) nor tree with hybrid edges (hasHybEdge) $(node.hasHybEdge), return the node.edge in order, unless a leaf is attached, then the edge attached to leaf is last");
        edge1 = nothing
        edge2 = nothing
        edge3 = nothing
        leaffound = false
        ind = 1
        for i in 1:3
            if(getOtherNode(node.edge[i],node).leaf)
                leaffound = true
                edge3 = node.edge[i]
                ind = i
                break
            end
        end
        if(leaffound)
            if(ind == 1)
                return node.edge[2], node.edge[3], edge3
            elseif(ind == 2)
                return node.edge[1], node.edge[3], edge3
            elseif(ind == 3)
                return node.edge[1], node.edge[2], edge3
            end
        else
            return node.edge[1], node.edge[2], node.edge[3]
        end
    end
end

# function to get the other two edges of a node
# besides the one specified
# it is called hybridEdges, but it not restricted to hybrid
function hybridEdges(node::Node, edge::Edge)
    size(node.edge,1) == 3 || error("node $(node.number) has $(size(node.edge,1)) edges instead of 3")
    edge1 = nothing
    edge2 = nothing
    for e in node.edge
        if(!isequal(e,edge))
            isa(edge1,Void) ? edge1 = e : edge2 = e
        end
    end
    return edge1,edge2
end

# function to remove an edge from a node
# warning: deletion is final, you can only
#          have edge back by pushing it again
# warning: if the edge removed is hybrid and node is tree,
#          node.hasHybEdge is set to false
#          assuming any tree node can only have one
#          one hybrid edge
function removeEdge!(node::Node,edge::Edge)
    index = 0
    try
        index = getIndexEdge(edge,node);
    catch e
        if isa(e, ErrorException)
            error("edge $(edge.number) not in node $(node.number)")
        end
    end
    deleteat!(node.edge,index);
    all((e->!e.hybrid), node.edge) ? node.hasHybEdge = false : node.hasHybEdge = true;
end

# function to remove a node from a edge
# warning: deletion is final, you can only
#          have node back by pushing it again
# warning: only removes node from edge, edge might still
#          be in node.edge
function removeNode!(node::Node,edge::Edge)
    index = 0
    try
        index = getIndexNode(edge,node);
    catch e
        if isa(e, ErrorException)
            error("node $(node.number) not in edge or strange edge with more than 2 nodes")
        end
    end
    deleteat!(edge.node,index);
end

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

# setLength
# warning: allows to change edge length for istIdentifiable=false
#          but issues a warning
# negative=true means it allows negative branch lengths (useful in qnet typeHyb=4)
function setLength!(edge::Edge, new_length::Number, negative::Bool)
    (negative || new_length >= 0) || error("length has to be nonnegative: $(new_length), cannot set to edge $(edge.number)")
    new_length >= -0.4054651081081644 || error("length can be negative, but not too negative (greater than -log(1.5)) or majorCF<0: new length is $(new_length)")
    #println("setting length $(new_length) to edge $(edge.number)")
    if(new_length > 10.0)
        new_length = 10.0;
    end
    edge.length = new_length;
    edge.y = exp(-new_length);
    edge.z = 1.0 - edge.y;
    #edge.istIdentifiable || warn("set edge length for edge $(edge.number) that is not identifiable")
    return nothing
end

"""
`setLength!(Edge,new length)`

set a new length for an object Edge. The new length needs to be positive.
For example, if you have a HybridNetwork object net, and do printEdges(net), you can see the list of Edges and their lengths. You can then change the length of the 3rd edge with setLength!(net.edge[3],1.2).
If `new length` is above 10, the value 10 will be used, as an upper limit
to coalescent units that can be reliably estimated.
"""
setLength!(edge::Edge, new_length::Number) = setLength!(edge, new_length, false)

"""
`setBranchLength!(Edge,new length)`

sets the length of an Edge object. The new length needs to be non-negative, or -1.0 to be interpreted as missing.
Example: if net is a HybridNetwork object, do printEdges(net) to see the list of all edges with their lengths. The length of the 3rd edge can be changed to 1.2 with setBranchLength!(net.edge[3],1.2). It can also be set to missing with setBranchLength!(net.edge[3],-1.0)
"""
function setBranchLength!(edge::Edge, new_length::Number)
    (new_length >= 0 || new_length == -1.0) || error("length $(new_length) has to be nonnegative or -1.0 (for missing).")
    edge.length = new_length;
    edge.y = exp(-new_length);
    edge.z = 1.0 - edge.y;
end

"""
`setGamma!(Edge,new gamma)`

Set inheritance probability γ for an edge, which must be a hybrid edge.
The new γ needs to be in [0,1]. The γ of the "partner" hybrid edge is changed
accordingly, to 1-γ. The field `isMajor` is also changed accordingly.
If the new γ is approximately 0.5, `Edge` is set to the major parent,
its partner is set to the minor parent.

If `net` is a HybridNetwork object, `printEdges(net)` will show the list of edges
and their γ's. The γ of the third hybrid edge (say) can be changed to 0.2 with
`setGamma!(net.edge[3],0.2)`.
This will automatically set γ of the partner hybrid edge to 0.8.
"""
# warning in the bad diamond/triangle cases because gamma is not identifiable
# changeOther = true, looks for the other hybrid edge and changes gamma too

setGamma!(edge::Edge, new_gamma::Float64) = setGamma!(edge, new_gamma, true)

function setGamma!(edge::Edge, new_gamma::Float64, changeOther::Bool)
    global DEBUG
    new_gamma >= 0.0 || error("gamma has to be positive: $(new_gamma)")
    new_gamma <= 1.0 || error("gamma has to be less than 1: $(new_gamma)")
    edge.hybrid || error("cannot change gamma in a tree edge");
    edge.isChild1 ? ind = 1 : ind = 2 ; # hybrid edge pointing at node 1 or 2
    node = edge.node[ind] # child of hybrid edge
    node.hybrid || warn("hybrid edge $(edge.number) not pointing at hybrid node")
    if(DEBUG)
        !node.isBadDiamondI || warn("bad diamond situation: gamma not identifiable")
    end
    edges = hybridEdges(node,edge)
    length(edges) == 2 || error("strange here: node $(node.number) should have 3 edges and it has $(length(edges)+1).")
    if(edges[1].hybrid && !edges[2].hybrid)
        ind = 1
    elseif(edges[2].hybrid && !edges[1].hybrid)
        ind = 2
    else
        error("strange hybrid node $(node.number) with only one hybrid edge or with three hybrid edges")
    end
    if(changeOther)
        if(!approxEq(new_gamma,0.5))
            edge.gamma = new_gamma;
            edge.isMajor = (new_gamma>0.5)
            edges[ind].gamma = 1.0 - new_gamma;
            edges[ind].isMajor = !edge.isMajor
        else #new gamma is 0.5 approximately
            edge.gamma = new_gamma
            edge.isMajor = true
            edges[ind].gamma = 1.0 - new_gamma
            edges[ind].isMajor = false
        end
    else
        if(!approxEq(new_gamma,0.5))
            edge.gamma = new_gamma;
            edge.isMajor = (new_gamma>0.5)
        else #new gamma is 0.5
            edge.gamma = new_gamma
            edge.isMajor = !edges[ind].isMajor
        end
    end
    return nothing
end

"""
`setGammaBLfromGammaz!(node, network)`

Update the γ values of the two sister hybrid edges in a bad diamond I, given the `gammaz` values
of their parent nodes, and updates the branch lengths t1 and t2 of their parent edges
(those across from the hybrid nodes), in such a way that t1=t2 and that these branch lengths
and γ values are consistent with the `gammaz` values in the network.
"""
# similar to the beginning of undoGammaz!, but does not update anything else than γ and t's.
# unlike undoGammaz!, does nothing (no error) if non-hybrid or not at bad diamond I.
function setGammaBLfromGammaz!(node::Node, net::HybridNetwork)
    if !node.isBadDiamondI || !node.hybrid
        return nothing
    end
    edge_maj, edge_min, tree_edge2 = hybridEdges(node);
    other_maj = getOtherNode(edge_maj,node);
    other_min = getOtherNode(edge_min,node);
    edgebla,tree_edge_incycle1,tree_edge = hybridEdges(other_min);
    edgebla,tree_edge_incycle2,tree_edge = hybridEdges(other_maj);
    if(approxEq(other_maj.gammaz,0.0) && approxEq(other_min.gammaz,0.0))
        edge_maj.gamma = 1.0 # γ and t could be anything if both gammaz are 0
        edge_min.gamma = 0.0 # will set t's to 0 and minor γ to 0.
        newt = 0.0
    else
        ((approxEq(other_min.gammaz,0.0) || other_min.gammaz >= 0.0) &&
         (approxEq(other_maj.gammaz,0.0) || other_maj.gammaz >= 0.0)    ) ||
            error("bad diamond I in node $(node.number) but missing (or <0) gammaz")
        ztotal = other_maj.gammaz + other_min.gammaz
        edge_maj.gamma = other_maj.gammaz / ztotal
        edge_min.gamma = other_min.gammaz / ztotal
        newt = -log(1-ztotal)
    end
    setLength!(tree_edge_incycle1,newt)
    setLength!(tree_edge_incycle2,newt)
end

function numTreeEdges(net::HybridNetwork)
    2*net.numTaxa - 3 + net.numHybrids
end

function numIntTreeEdges(net::HybridNetwork)
    2*net.numTaxa - 3 + net.numHybrids - net.numTaxa
end

# function to get the partition where an edge is
# returns the index of the partition, or error if not found
# better to return the index than the partition itself, because we need the index
# to use splice and delete it from net.partition later on
# cycle: is the number to look for partition on that cycle only
function whichPartition(net::HybridNetwork,edge::Edge,cycle::Integer)
    global DEBUG
    !edge.hybrid || error("edge $(edge.number) is hybrid so it cannot be in any partition")
    edge.inCycle == -1 || error("edge $(edge.number) is in cycle $(edge.inCycle) so it cannot be in any partition")
    DEBUG && println("search partition for edge $(edge.number) in cycle $(cycle)")
    in(edge,net.edge) || error("edge $(edge.number) is not in net.edge")
    for i in 1:length(net.partition)
        DEBUG && println("looking for edge $(edge.number) in partition $(i): $([e.number for e in net.partition[i].edges])")
        if(in(cycle,net.partition[i].cycle))
            DEBUG && println("looking for edge $(edge.number) in partition $(i), with cycle $(cycle): $([e.number for e in net.partition[i].edges])")
            if(in(edge,net.partition[i].edges))
                DEBUG && println("partition for edge $(edge.number) is $([e.number for e in net.partition[i].edges])")
                return i
            end
        end
    end
    DEBUG && printPartitions(net)
    error("edge $(edge.number) is not hybrid, nor part of any cycle, and it is not in any partition")
end

# function to get the partition where an edge is
# returns the index of the partition, or error if not found
# better to return the index than the partition itself, because we need the index
# to use splice and delete it from net.partition later on
function whichPartition(net::HybridNetwork,edge::Edge)
    global DEBUG
    !edge.hybrid || error("edge $(edge.number) is hybrid so it cannot be in any partition")
    edge.inCycle == -1 || error("edge $(edge.number) is in cycle $(edge.inCycle) so it cannot be in any partition")
    DEBUG && println("search partition for edge $(edge.number) without knowing its cycle")
    in(edge,net.edge) || error("edge $(edge.number) is not in net.edge")
    for i in 1:length(net.partition)
        DEBUG && println("looking for edge $(edge.number) in partition $(i): $([e.number for e in net.partition[i].edges])")
        if(in(edge,net.partition[i].edges))
            DEBUG && println("partition for edge $(edge.number) is $([e.number for e in net.partition[i].edges])")
            return i
        end
    end
    DEBUG && printPartitions(net)
    error("edge $(edge.number) is not hybrid, nor part of any cycle, and it is not in any partition")
end

# function that will print the partition of net
function printPartitions(net::HybridNetwork)
    println("partition.cycle\t partition.edges")
    for p in net.partition
        println("$(p.cycle)\t\t $([e.number for e in p.edges])")
    end
end

# function to find if a given partition is in net.partition
function isPartitionInNet(net::HybridNetwork,desc::Vector{Edge},cycle::Vector{Int})
    for p in net.partition
        if(sort(cycle) == sort(p.cycle))
            if(sort([e.number for e in desc]) == sort([e.number for e in p.edges]))
                return true
            end
        end
    end
    return false
end

# function to check that everything matches in a network
# in particular, cycles, partitions and containRoot
# fixit: need to add check on identification of bad diamonds, triangles
# and correct computation of gammaz
# light=true: it will not collapse with nodes with 2 edges, will return a flag of true
# returns true if found egde with BL -1.0 (only when light=true, ow error)
# added checkPartition for undirectedOtherNetworks that do not need correct hybrid node number
function checkNet(net::HybridNetwork, light::Bool; checkPartition=true::Bool)
    global DEBUG
    DEBUG && println("checking net")
    net.numHybrids == length(net.hybrid) || error("discrepant number on net.numHybrids (net.numHybrids) and net.hybrid length $(length(net.hybrid))")
    net.numTaxa == length(net.leaf) || error("discrepant number on net.numTaxa (net.numTaxa) and net.leaf length $(length(net.leaf))")
    net.numNodes == length(net.node) || error("discrepant number on net.numNodes (net.numNodes) and net.node length $(length(net.node))")
    net.numEdges == length(net.edge) || error("discrepant number on net.numEdges (net.numEdges) and net.edge length $(length(net.edge))")
    if(isTree(net))
        all(x->x.containRoot,net.edge) || error("net is a tree, but not all edges can contain root")
        all(x->x.isMajor,net.edge) || error("net is a tree, but not all edges are major")
        all(x->!(x.hybrid),net.edge) || error("net is a tree, but not all edges are tree")
        all(x->!(x.hybrid),net.node) || error("net is a tree, but not all nodes are tree")
        all(x->!(x.hasHybEdge),net.node) || error("net is a tree, but not all nodes hasHybEdge=false")
        all(x->(x.gamma == 1.0 ? true : false),net.edge) || error("net is a tree, but not all edges have gamma 1.0")
    end
    for h in net.hybrid
        if(isBadTriangle(h))
            DEBUG && println("hybrid $(h.number) is very bad triangle")
            net.hasVeryBadTriangle || error("hybrid node $(h.number) is very bad triangle, but net.hasVeryBadTriangle is $(net.hasVeryBadTriangle)")
            h.isVeryBadTriangle || h.isExtBadTriangle || error("hybrid node $(h.number) is very bad triangle but it does not know it")
        end
        nocycle,edges,nodes = identifyInCycle(net,h)
        for e in edges
            e.inCycle == h.number || error("edge $(e.number) is in cycle of hybrid node $(h.number) but its inCycle attribute is $(e.inCycle)")
            if(e.length == -1.0)
                if(light)
                    return true
                else
                    error("found edge with BL -1.0")
                end
            end
            if(e.hybrid)
                !e.containRoot || error("hybrid edge $(e.number) should not contain root") # fixit: disagree
                o = getOtherNode(e,h)
                o.hasHybEdge || error("found node $(o.number) attached to hybrid edge but hasHybEdge=$(o.hasHybEdge)")
            end
        end
        for n in nodes
            n.inCycle == h.number || error("node $(n.number) is in cycle of hybrid node $(h.number) but its inCycle attribute is $(n.inCycle)")
            e1,e2,e3 = hybridEdges(n)
            i = 0
            for e in [e1,e2,e3]
                if(isa(e,Void) && h.k != 2)
                    error("edge found that is Void, and hybrid node $(h.number) k is $(h.k). edge as nothing can only happen when k=2")
                elseif(!isa(e,Void))
                    if(e.inCycle == -1)
                        i += 1
                        desc = [e]
                        cycleNum = [h.number]
                        getDescendants!(getOtherNode(e,n),e,desc,cycleNum)
                        if(checkPartition && !isPartitionInNet(net,desc,cycleNum))
                            printPartitions(net)
                            error("partition with cycle $(cycleNum) and edges $([e.number for e in desc]) not found in net.partition")
                        end
                    end
                end
            end
            i == 1 || error("strange node $(n.number) incycle $(h.number) but with $(i) edges not in cycle, should be only one")
            edgesRoot = identifyContainRoot(net,h)
            for edge in edgesRoot
                if(edge.containRoot)
                    DEBUG && printEverything(net)
                    error("edge $(edge.number) should not contain root")
                end
            end
        end
    end
    for n in net.node
        if(n.leaf)
            length(n.edge) == 1 || error("leaf $(n.number) with $(length(n.edge)) edges instead of 1")
        else
            if(light)
                if(length(n.edge) != 3)
                    DEBUG && warn("node $(n.number) with $(length(n.edge)) edges instead of 3")
                    return true
                end
            else
                length(n.edge) == 3 || error("node $(n.number) with $(length(n.edge)) edges instead of 3")
            end
        end
    end
    for e in net.edge
        if(e.length == -1.0)
            if(light)
                return true
            else
                error("edge found with BL -1.0")
            end
        end
    end
    DEBUG && println("no errors in checking net")
    return false
end

checkNet(net::HybridNetwork) = checkNet(net, false)

# function to print everything for a given net
# this is used a lot inside snaq when DEBUG is true, so need to use level1 attributes
# and not change the network: with writeTopologyLevel1
function printEverything(net::HybridNetwork)
    printEdges(net)
    printNodes(net)
    printPartitions(net)
    println("$(writeTopologyLevel1(net))")
end

# function to check if a node is very or ext bad triangle
function isBadTriangle(node::Node)
    node.hybrid || error("cannot check if node $(node.number) is very bad triangle because it is not hybrid")
    if(node.k == 3)
        edgemaj, edgemin, treeedge = hybridEdges(node)
        othermaj = getOtherNode(edgemaj,node)
        othermin = getOtherNode(edgemin,node)
        treenode = getOtherNode(treeedge,node)
        edges1 = hybridEdges(othermaj)
        o1 = getOtherNode(edges1[3],othermaj)
        edges2 = hybridEdges(othermin)
        o2 = getOtherNode(edges2[3],othermin)
        leaves = sum([n.leaf ? 1 : 0 for n in [treenode,o1,o2]])
        if(leaves == 1 || leaves == 2)
            return true
        else
            return false
        end
    else
        return false
    end
end

# function to check if a partition is already in net.partition
# used in updatePartition
function isPartitionInNet(net::HybridNetwork,partition::Partition)
    if(isempty(net.partition))
        return false
    end
    for p in net.partition
        cycle = isempty(setdiff(p.cycle,partition.cycle)) && isempty(setdiff(partition.cycle,p.cycle))
        edges = isempty(setdiff([n.number for n in p.edges],[n.number for n in partition.edges])) && isempty(setdiff([n.number for n in partition.edges],[n.number for n in p.edges]))
        if(cycle && edges)
            return true
        end
    end
    return false
end

# function to switch a hybrid node in a network to another node in the cycle
function switchHybridNode!(net::HybridNetwork, hybrid::Node, newHybrid::Node)
    hybrid.hybrid || error("node $(hybrid.number) has to be hybrid to switch to a different hybrid")
    newHybrid.inCycle == hybrid.number || error("new hybrid needs to be in the cycle of old hybrid: $(hybrid.number)")
    !newHybrid.hybrid || error("strange hybrid node $(newHybrid.number) in cycle of another hybrid $(hybrid.number)")
    newHybrid.hybrid = true
    newHybrid.hasHybEdge = true
    newHybrid.name = hybrid.name
    pushHybrid!(net,newHybrid)
    makeNodeTree!(net,hybrid)
end

"""
`assignhybridnames!(net)`

Assign names to hybrid nodes in the network `net`. Hybrid nodes with an empty `name` field ("")
are modified with a name that does not conflict with other hybrid names in the network.
The preferred name is "#H3" if the node number is 3 or -3, but an index other than 3 would be
used if "#H3" were the name of another hybrid node already.

If two hybrid nodes have non-empty and equal names, the name of one of them is changed and
re-assigned as described above (with a warning).
"""
function assignhybridnames!(net::HybridNetwork)
    hybnum = Int[]  # indices 'i' in hybrid names: #Hi
    # first: go through *all* existing non-empty names
    for ih in 1:length(net.hybrid)
        lab = net.hybrid[ih].name
        lab != "" || continue # do nothing if label is missing
        jh = findfirst([net.hybrid[j].name for j in 1:ih-1], lab)
        if jh > 0 # set repeated names to ""
            warn("hybrid nodes $(net.hybrid[ih].number) and $(net.hybrid[jh].number) have the same label: $lab. Will change the name of the former.")
            net.hybrid[ih].name = ""
        else
            m = match(r"^#H(\d+)$", lab)
            if m != nothing # make full list of existing indices "i" in #Hi
                push!(hybnum, parse(Int, m[1]))
            end
        end
    end
    # second: assign empty names to #Hi
    hnext = 1
    for ih in 1:length(net.hybrid)
        net.hybrid[ih].name == "" || continue # do nothing if non-empty label
        hnum = abs(net.hybrid[ih].number)
        while in(hnum, hybnum)
            hnum = hnext  # not efficient, but rare
            hnext += 1    # and okay on small networks
        end
        push!(hybnum, hnum)
        net.hybrid[ih].name = "#H$hnum"
    end
end

"""
    sorttaxa!(DataFrame, columns)

Reorder the 4 taxa and reorders the observed concordance factors accordingly, on each row of
the data frame. If `columns` is committed, taxon names are assumed to be in columns 1-4 and
CFs are assumed to be in columns 5-6 with quartets in this order: 12_34, 13_24, 14_23.
Does **not** reorder credibility interval values, if present.

    sorttaxa!(DataCF)
    sorttaxa!(Quartet, permutation_tax, permutation_cf)

Reorder the 4 taxa in each element of the DataCF `quartet`. For a given Quartet,
reorder the 4 taxa in its fields `taxon` and `qnet.quartetTaxon` (if non-empty)
and reorder the 3 concordance values accordingly, in `obsCF` and `qnet.expCF`.

`permutation_tax` and `permutation_cf` should be vectors of short integers (Int8) of length 4 and 3
respectively, whose memory allocation gets reused. Their length is *not checked*.
"""
# qnet.names unchanged: order of taxon names here relates to the order of nodes in the network

function sorttaxa!(dat::DataCF)
    ptax = Array{Int8}(4) # to hold the sort permutations
    pCF  = Array{Int8}(3)
    for q in dat.quartet
        sorttaxa!(q, ptax, pCF)
    end
end

function sorttaxa!(df::DataFrame, co=Int[]::Vector{Int})
    if length(co)==0
        co = collect(1:7)
    end
    length(co) > 6 || error("column vector must be of length 7 or more")
    ptax = Array{Int8}(4)
    pCF  = Array{Int8}(3)
    taxnam = Array{eltype(df[co[1]])}(4)
    for i in 1:size(df,1)
        for j=1:4 taxnam[j] = df[i,co[j]]; end
        sortperm!(ptax, taxnam)
        sorttaxaCFperm!(pCF, ptax) # update permutation pCF according to taxon permutation
        df[i,co[1]], df[i,co[2]], df[i,co[3]], df[i,co[4]] = taxnam[ptax[1]], taxnam[ptax[2]], taxnam[ptax[3]], taxnam[ptax[4]]
        df[i,co[5]], df[i,co[6]], df[i,co[7]] = df[i,co[pCF[1]+4]], df[i,co[pCF[2]+4]], df[i,co[pCF[3]+4]]
    end
    return df
end

function sorttaxa!(qua::Quartet, ptax::Vector{Int8}, pCF::Vector{Int8})
    qt = qua.taxon
    if length(qt)==4
        sortperm!(ptax, qt)
        sorttaxaCFperm!(pCF, ptax) # update permutation pCF accordingly
        qt[1], qt[2], qt[3], qt[4] = qt[ptax[1]], qt[ptax[2]], qt[ptax[3]], qt[ptax[4]]
        qua.obsCF[1], qua.obsCF[2], qua.obsCF[3] = qua.obsCF[pCF[1]], qua.obsCF[pCF[2]], qua.obsCF[pCF[3]]
        # do *NOT* modify qua.qnet.quartetTaxon: it points to the same array as qua.taxon
        eCF = qua.qnet.expCF
        if length(eCF)==3
            eCF[1], eCF[2], eCF[3] = eCF[pCF[1]], eCF[pCF[2]], eCF[pCF[3]]
        end
    elseif length(qt)!=0
        error("Quartet with $(length(qt)) taxa")
    end
    return qua
end

# find permutation pCF of the 3 CF values: 12_34, 13_24, 14_23. 3!=6 possible permutations
# ptax = one of 4!=24 possible permutations on the 4 taxon names
# kernel: pCF = identity if ptax = 1234, 2143, 3412 or 4321
# very long code, but to minimize equality checks at run time
function sorttaxaCFperm!(pcf::Vector{Int8}, ptax::Vector{Int8})
    if ptax[1]==1
        if     ptax[2]==2
            pcf[1]=1
            if  ptax[3]==3 # ptax = 1,2,3,4
                pcf[2]=2; pcf[3]=3
            else           # ptax = 1,2,4,3
                pcf[2]=3; pcf[3]=2
            end
        elseif ptax[2]==3
            pcf[1]=2
            if  ptax[3]==2 # ptax = 1,3,2,4
                pcf[2]=1; pcf[3]=3
            else           # ptax = 1,3,4,2
                pcf[2]=3; pcf[3]=1
            end
        else # ptax[2]==4
            pcf[1]=3
            if  ptax[3]==2 # ptax = 1,4,2,3
                pcf[2]=1; pcf[3]=2
            else           # ptax = 1,4,3,2
                pcf[2]=2; pcf[3]=1
            end
        end
    elseif ptax[1]==2
        if     ptax[2]==1
            pcf[1]=1
            if  ptax[3]==4 # ptax = 2,1,4,3
                pcf[2]=2; pcf[3]=3
            else           # ptax = 2,1,3,4
                pcf[2]=3; pcf[3]=2
            end
        elseif ptax[2]==4
            pcf[1]=2
            if  ptax[3]==1 # ptax = 2,4,1,3
                pcf[2]=1; pcf[3]=3
            else           # ptax = 2,4,3,1
                pcf[2]=3; pcf[3]=1
            end
        else # ptax[2]==3
            pcf[1]=3
            if  ptax[3]==1 # ptax = 2,3,1,4
                pcf[2]=1; pcf[3]=2
            else           # ptax = 2,3,4,1
                pcf[2]=2; pcf[3]=1
            end
        end
    elseif ptax[1]==3
        if     ptax[2]==4
            pcf[1]=1
            if  ptax[3]==1 # ptax = 3,4,1,2
                pcf[2]=2; pcf[3]=3
            else           # ptax = 3,4,2,1
                pcf[2]=3; pcf[3]=2
            end
        elseif ptax[2]==1
            pcf[1]=2
            if  ptax[3]==4 # ptax = 3,1,4,2
                pcf[2]=1; pcf[3]=3
            else           # ptax = 3,1,2,4
                pcf[2]=3; pcf[3]=1
            end
        else # ptax[2]==2
            pcf[1]=3
            if  ptax[3]==4 # ptax = 3,2,4,1
                pcf[2]=1; pcf[3]=2
            else           # ptax = 3,2,1,4
                pcf[2]=2; pcf[3]=1
            end
        end
    else # ptax[1]==4
        if     ptax[2]==3
            pcf[1]=1
            if  ptax[3]==2 # ptax = 4,3,2,1
                pcf[2]=2; pcf[3]=3
            else           # ptax = 4,3,1,2
                pcf[2]=3; pcf[3]=2
            end
        elseif ptax[2]==2
            pcf[1]=2
            if  ptax[3]==3 # ptax = 4,2,3,1
                pcf[2]=1; pcf[3]=3
            else           # ptax = 4,2,1,3
                pcf[2]=3; pcf[3]=1
            end
        else # ptax[2]==1
            pcf[1]=3
            if  ptax[3]==3 # ptax = 4,1,3,2
                pcf[2]=1; pcf[3]=2
            else           # ptax = 4,1,2,3
                pcf[2]=2; pcf[3]=1
            end
        end
    end
end


text:
"""
`module  JAC.PhotoIonizationAutoion`
    ... a submodel of JAC that contains all methods for computing photo-excitation-autoionization cross sections and rates.
"""
module PhotoIonizationAutoion

    using ..AutoIonization, ..Basics, ..ManyElectron, ..Radial, ..PhotoEmission, ..TableStrings

    """
    `struct  PhotoIonizationAutoion.Settings`
        ... defines a type for the details and parameters of computing photon-impact excitation-autoionization pathways
            |i(N)>  --> |m(N)>  --> |f(N-1)>.

        + multipoles              ::Array{EmMultipole,1}               ... Specifies the multipoles of the radiation field that are to be included.
        + gauges                  ::Array{UseGauge,1}                  ... Specifies the gauges to be included into the computations.
        + printBefore             ::Bool                               ... True, if all energies and lines are printed before their evaluation.
        + selectPathways          ::Bool                               ... True if particular pathways are selected for the computations.
        + incidentStokes          ::ExpStokes               ... Stokes parameters of the incident radiation.
        + selectedPathways        ::Array{Tuple{Int64,Int64,Int64},1}  ... List of list of pathways, given by tupels (inital, inmediate, final).
        + maxKappa                ::Int64                              ... Maximum kappa value of partial waves to be included.
    """
    struct Settings
        multipoles                ::Array{EmMultipole,1}
        gauges                    ::Array{UseGauge,1}
        printBefore    ::Bool
        selectPathways            ::Bool
        selectedPathways          ::Array{Tuple{Int64,Int64,Int64},1}
        maxKappa                  ::Int64
    end

    """
    `PhotoIonizationAutoion.Settings()`  ... constructor for the default values of photon-impact excitation-autoionizaton settings.
    """
    function Settings()
        Settings( Basics.EmMultipole[], Basics.UseGauge[], false,  false, Tuple{Int64,Int64,Int64}[], 0)
    end

    # `Base.show(io::IO, settings::PhotoIonizationAutoion.Settings)`
    #		 ... prepares a proper printout of the variable settings::PhotoIonizationAutoion.Settings.
    function Base.show(io::IO, settings::PhotoIonizationAutoion.Settings)
        println(io, "multipoles:              $(settings.multipoles)  ")
        println(io, "gauges:                  $(settings.gauges)  ")
        println(io, "printBefore:  $(settings.printBefore)  ")
        println(io, "selectPathways:          $(settings.selectPathways)  ")
        println(io, "selectedPathways:        $(settings.selectedPathways)  ")
        println(io, "maxKappa:                $(settings.maxKappa)  ")
    end

    """
    `struct  PhotoIonizationAutoion.Channel`
        ... defines a type for a photon-impact excitaton & autoionization channel that specifies all quantum numbers, phases and
            amplitudes.

        + excitationChannel  ::PhotoEmission.Channel       ... Channel that describes the photon-impact excitation process.
        + augerChannel       ::AutoIonization.Channel      ... Channel that describes the subsequent Auger/autoionization process.
    """
    struct  Channel
        excitationChannel    ::PhotoEmission.Channel
        augerChannel         ::AutoIonization.Channel
    end

    """
    `struct  PhotoIonizationAutoion.Pathway`  ... defines a type for a photon-impact excitation pathway that may include the definition
                                                      of different excitation and autoionization channels and their corresponding amplitudes.

        + initialLevel        ::Level           ... initial-(state) level
        + intermediateLevel   ::Level           ... intermediate-(state) level
        + finalLevel          ::Level           ... final-(state) level
        + photonEnergy        ::Float64         ... energy of the (incoming) electron
        + electronEnergy      ::Float64         ... energy of the (finally outgoing, scattered) electron
        + crossSection        ::EmProperty      ... total cross section of this pathway
        + hasChannels         ::Bool            ... Determines whether the individual excitation and autoionization channels
                                                    are defined in terms of their multipole, gauge, free-electron kappa,
                                                    phases and the total angular momentum/parity as well as the amplitude, or not.
        + channels            ::Array{PhotoIonizationAutoion.Channel,1}     ... List of channels of this pathway.
    """
    struct  Pathway
        initialLevel          ::Level
        intermediateLevel     ::Level
        finalLevel            ::Level
        photonEnergy          ::Float64
        electronEnergy        ::Float64
        crossSection          ::EmProperty
        hasChannels           ::Bool
        channels              ::Array{PhotoIonizationAutoion.Channel,1}
    end

    """
    `PhotoIonizationAutoion.Pathway()`
        ... 'empty' constructor for an photon-impact excitation-autoionization pathway between a specified initial, intermediate and
            final level.
    """
    function Pathway()
        Pathway(Level(), Level(), Level(), 0., 0., 0., false, PhotoIonizationAutoion.Channel[] )
    end

    # `Base.show(io::IO, pathway::PhotoIonizationAutoion.Pathway)`
    #		 ... prepares a proper printout of the variable pathway::PhotoIonizationAutoion.Pathway.
    function Base.show(io::IO, pathway::PhotoIonizationAutoion.Pathway)
        println(io, "initialLevel:               $(pathway.initialLevel)  ")
        println(io, "intermediateLevel:          $(pathway.intermediateLevel)  ")
        println(io, "finalLevel:                 $(pathway.finalLevel)  ")
        println(io, "photonEnergy                $(pathway.photonEnergy)  ")
        println(io, "electronEnergy              $(pathway.electronEnergy)  ")
        println(io, "crossSection:               $(pathway.crossSection)  ")
        println(io, "hasChannels:                $(pathway.hasChannels)  ")
        println(io, "channels:                   $(pathway.channels)  ")
    end

    """
    `PhotoIonizationAutoion.computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet,
                                            grid::Radial.Grid, settings::PhotoIonizationAutoion.Settings; output=true)`
        ... to compute the photo-excitation-fluorescence amplitudes and all properties as requested by the given settings. A list of
            lines::Array{PhotoIonizationAutoion.Lines} is returned.
    """
    function  computePathways(finalMultiplet::Multiplet, intermediateMultiplet::Multiplet, initialMultiplet::Multiplet, grid::Radial.Grid,
                              settings::PhotoIonizationAutoion.Settings; output=true)
        println("")
        printstyled("PhotoIonizationAutoion.computePathways(): The computation of photo-excitation-fluorescence amplitudes starts now ... \n", color=:light_green)
        printstyled("-------------------------------------------------------------------------------------------------------------------- \n", color=:light_green)
        println("")
        #
        println("Not yet implemented: Data structures and properties still need to be worked out.")
        #
        pathways = "Not yet implemented !"
        #
        if    output    return( pathways )
        else            return( nothing )
        end
    end

end # module

ext: jl


text: # File empty on purpose. Tests will be added shortly.

max_stars_repo_name: PCwarlord/Julia-GPU-CPU-balancer

max_issues_repo_name: PCwarlord/Julia-GPU-CPU-balancer

max_forks_repo_name: PCwarlord/Julia-GPU-CPU-balancer


text:
const SFACT = Symbol("@fact")
const SFACT_THROWS = Symbol("@fact_throws")

const COMPARISON_DICT = Dict(
:exactly               => Symbol("==="),
:less_than             => Symbol("<"),
:less_than_or_equal    => Symbol("<="),
:greater_than          => Symbol(">"),
:greater_than_or_equal => Symbol(">="),
)

const PROPERTIES = [:isnan, :isinf, :isfinite, :isodd, :iseven, :isempty]

"""
    fc2ts(x)

Convert expression like object x from FactCheck to TestSet.
"""
fc2ts(x) = x
fc2ts(x::Expr) = fc2ts(x, extrait(x))

function fc2ts(ex::Expr, ::MacroCall{SFACT})
     ex |> fact2testcore |> testex
end

function fc2ts(ex::Expr, ::Any)
    args = map(fc2ts, ex.args)
    return Expr(ex.head, args...)
end

function fc2ts(ex::Expr, ::MacroCall{SFACT_THROWS})
    args = ex.args[2:end]
    if length(args) == 1
        except = :Exception
        code = args[1]
    else
        @assert length(args) == 2
        except, code = args
    end
    :(@test_throws $except $code)
end

function fc2ts(ex::Expr, ::Call{:context})
    name, arg, body = name_arg_body(ex)
    body = fc2ts(body)
    :(@testset $arg $body)
end

function fc2ts(ex::Expr, ::Call{:facts})
    name, arg, body = name_arg_body(ex)
    body = fc2ts(body)
    :(@testset $arg $body)
end

function name_arg_body(doblock::Expr)
    @match doblock begin
        (name_)(arg_) do
            body_
        end => (name, arg, body)
    end
end

fact2testcore(ex) = lhs_rhs2testcore(lhs_rhs(ex)...)

function lhs_rhs(ex)
    @match ex begin
        @fact(lhs_ --> rhs_) => (lhs, rhs)
    end
end

testex(args...) = Expr(:macrocall, Symbol("@test"), args...)

lhs_rhs2testcore(lhs, rhs) = lhs_rhs2testcore(lhs, rhs, extrait(rhs))
lhs_rhs2testcore(lhs, rhs, rtrait) = :($lhs == $rhs)
lhs_rhs2testcore(lhs, rhs::Bool, ::NoExpr) = rhs ? lhs : :(!$lhs)

function lhs_rhs2testcore(lhs, rhs::Symbol, ::NoExpr)
    if rhs in PROPERTIES
        return :($rhs($lhs))
    else
        return :($lhs == $rhs)
    end
end

function lhs_rhs2testcore(lhs, rhs::Expr, ::Call{:not})
    @assert rhs.args[1] == :not
    rhs_new = callarg_unique(rhs)
    ex = lhs_rhs2testcore(lhs, rhs_new)
    :(!$ex)
end

atol2kw(x) = Expr(:kw, :atol, x)
function atol2kw(x::Expr)
    x.head == :kw && return x
    return Expr(:kw, :atol, x)
end

function lhs_rhs2testcore(lhs, rhs::Expr, ::Call{:roughly})
    args = callargs(rhs)
    arg = args[1]
    kw = map(atol2kw, args[2:end])
    if length(args) == 1
        return :($lhs ≈ $arg)
    else
        return Expr(:call, :isapprox, lhs, arg, kw...)
    end
end

function lhs_rhs2testcore{f}(lhs, rhs::Expr, ::Call{f})
    if f in keys(COMPARISON_DICT)
        s = COMPARISON_DICT[f]
        rhs_new = callarg_unique(rhs)
        return Expr(:call, s, lhs, rhs_new)
    else
        return :($lhs == $rhs)
    end
end

function lhs_rhs2testcore(lhs, rhs::Expr, ::Lambda)
    Expr(:call, rhs, lhs)
end


text:
"""
    Observable

Represents a physical quantity that can be observed, which is
calculated as the matrix element of an operator between two
configurations. All physical observables are real.
"""
mutable struct Observable{T,B<:Basis,Equations<:AbstractVector{<:AtomicOrbitalEquation}}
    equations::Equations
    tmp::RadialOrbital{T,B}
end

"""
    Observable(operator, atom, overlaps, integrals, integral_map, selector
               [; double_counted=false])

Construct an observable corresponding the `operator` acting on `atom`;
if `double_counted`, only return those terms that would be
double-counted, otherwise return the normal observable
equations. `overlaps` is a list of non-orthogonal, `integrals` a list
of common integrals, and `integral_map` is a mapping from symbolic
integrals to [`OrbitalIntegral`](@ref)s. `selector` selects those
orbitals not modelled by the potential, e.g. all orbitals in case of a
nuclear potential, fewer in case of a
pseudopotential. `double_counted` can be used to only return those
terms of the sum that would be double-counted if simply summing over
orbital contributions (applicable to two-body operators only). Half
the result of the double-counted term is then subtracted the sum over
orbital contributions; this makes the expression slightly more
symmetric in orbital space than if the sum is derived avoiding
double-counting.
"""
function Observable(operator::QuantumOperator, atom::A,
                    overlaps::Vector{<:OrbitalOverlap},
                    integrals::Vector,
                    integral_map::Dict{Any,Int},
                    symmetries::Dict,
                    selector::Function;
                    double_counted::Bool=false,
                    kwargs...) where {T,A<:Atom{T}}
    M = Matrix(operator, selector.(atom.configurations), overlaps)

    m,n = size(M)

    if double_counted
        numbodies(operator) > 2 &&
            throw(ArgumentError("Cannot construct double-counted observables for N-body operators with N > 2."))
        M *= one(T)/2
        # We need to filter out those terms that are due to zero- and
        # one-body interactions, since they would not be
        # double-counted.
        for i = 1:m
            for j = 1:n
                iszero(M[i,j]) && continue
                M[i,j] = NBodyMatrixElement(filter(t -> any(f -> numbodies(f) > 1,
                                                           t.factors),
                                                   M[i,j].terms))
            end
        end
    end

    symbolic_integrals = []

    equation_system = map(atom.orbitals) do orbital
        corb = conj(orbital)
        equation = orbital_equation(M, corb, symbolic_integrals)

        # We need to filter out those terms that are dependent on
        # corb, to avoid double-counting when computing the
        # observables.
        for i = 1:m
            for j = 1:n
                (iszero(M[i,j]) || double_counted) && continue
                M[i,j] = NBodyMatrixElement(filter(t -> !isdependent(t, corb),
                                                   M[i,j].terms))
            end
        end
        equation
    end

    # TODO: Should observables really project out symmetry orbitals?
    eqs = generate_atomic_orbital_equations(
        atom, MCEquationSystem(equation_system, symbolic_integrals),
        integrals, integral_map,
        symmetries; kwargs...)

    tmp = similar(first(eqs).ϕ)
    tmp.args[2] .= false # To clear any NaNs

    Observable(eqs, tmp)
end

realifreal(::Type{R}, v) where {R<:Real} = real(v)
realifreal(::Type{R}, v) where {R} = v

"""
    observe!(A::M, o::Observable)

Compute the observable `o` between all configurations and store the
results as matrix elements of `A`.
"""
function observe!(A::M, o::Observable{T}) where {R,M<:AbstractMatrix{R},T}
    A .= zero(T)
    for eq in o.equations
        for term in eq.hamiltonian.terms
            materialize!(MulAdd(coefficient(term), term.A, eq.ϕ, zero(T), o.tmp))
            A[term.i,term.j] += realifreal(R, (materialize(applied(*, eq.ϕ', o.tmp))))
        end
    end
    A
end

"""
    observe!(A::M, o::Observable, atom::Atom)

Compute the observable `o` between all configurations, store the
results as matrix elements of `A`, and finally contract with respect
to the mixing coefficients of `atom` and return the result as a
scalar.
"""
function observe!(A::M, o::Observable{T}, atom::Atom) where {R,M<:AbstractMatrix{R},T}
    observe!(A, o)
    atom.mix_coeffs'A*atom.mix_coeffs
end


text:
# These `deserialize` methods are used to be able to deserialize intervals using the
# structure that was used before Intervals 1.3.

function Serialization.deserialize(s::AbstractSerializer, ::Type{Interval{T}}) where T
    left = deserialize(s)
    right = deserialize(s)
    inclusivity = deserialize(s)

    L = bound_type(first(inclusivity))
    R = bound_type(last(inclusivity))

    return Interval{T,L,R}(left, right)
end

function Serialization.deserialize(s::AbstractSerializer, ::Type{AnchoredInterval{P,T}}) where {P,T}
    anchor = deserialize(s)
    inclusivity = deserialize(s)

    L = bound_type(first(inclusivity))
    R = bound_type(last(inclusivity))

    return AnchoredInterval{P,T,L,R}(anchor)
end

# Works around an issue where a `StepRangeLen` tries to convert the step to an
# `AbstractInterval`.
# Issue introduced in https://github.com/JuliaLang/julia/pull/23194 (a698230daf6)
# and fixed in https://github.com/JuliaLang/julia/pull/34412 (fca037a162)
if v"0.7.0-DEV.1399" <= VERSION < v"1.5.0-DEV.127"
    Base.step(r::StepRangeLen{<:AbstractInterval}) = r.step
end


text:
# LOGICCIRCUITS LIBRARY ROOT

module LogicCircuits

using Reexport

include("Utils/Utils.jl")
@reexport using .Utils

include("abstract_logic_nodes.jl")
include("bit_circuit.jl")

include("queries/queries.jl")
include("queries/satisfies.jl")
include("queries/satisfies_flow.jl")

include("plain_logic_nodes.jl")

include("structured/abstract_vtrees.jl")
include("structured/plain_vtrees.jl")
include("structured/structured_logic_nodes.jl")

include("transformations.jl")

include("bdd/bdds.jl")

include("sdd/sdds.jl")
include("sdd/sdd_functions.jl")
include("sdd/apply.jl")

include("io/io.jl")

end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
function calculate_sum(upper_limit)
    s = 0
    for i in 1 : upper_limit - 1
        if i % 3 == 0 || i % 5 == 0
            s += i
        end
    end
    println(s)
end

calculate_sum(1000)

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
function String{P}(x::ArbFloat{P}, ndigits::Int, flags::UInt)
   n = max(1,min(abs(ndigits), floor(Int, P*0.3010299956639811952137)))
   cstr = ccall(@libarb(arb_get_str), Ptr{UInt8}, (Ptr{ArbFloat}, Int, UInt), &x, n, flags)
   s = unsafe_string(cstr)
   ccall(@libflint(flint_free), Void, (Ptr{UInt8},), cstr)
   s
end

function String{P}(x::ArbFloat{P}, flags::UInt)
   n = floor(Int, P*0.3010299956639811952137)
   cstr = ccall(@libarb(arb_get_str), Ptr{UInt8}, (Ptr{ArbFloat}, Int, UInt), &x, n, flags)
   s = unsafe_string(cstr)
   ccall(@libflint(flint_free), Void, (Ptr{UInt8},), cstr)
   s
end

function string{P}(x::ArbFloat{P}, ndigits::Int)
   # n=trunc(abs(log(upperbound(x)-lowerbound(x))/log(2))) just the good bits
   s = String(x, ndigits, UInt(2)) # midpoint only (within 1ulp), RoundNearest
   s
end

function string{P}(x::ArbFloat{P})
   # n=trunc(abs(log(upperbound(x)-lowerbound(x))/log(2))) just the good bits
   s = String(x,UInt(2)) # midpoint only (within 1ulp), RoundNearest
   s
end

function stringTrimmed{P}(x::ArbFloat{P}, ndigitsremoved::Int)
   n = max(0, ndigitsremoved)
   n = max(1, floor(Int, P*0.3010299956639811952137) - n)
   cstr = ccall(@libarb(arb_get_str), Ptr{UInt8}, (Ptr{ArbFloat}, Int, UInt), &x, n, UInt(2))
   s = unsafe_string(cstr)
   ccall(@libflint(flint_free), Void, (Ptr{UInt8},), cstr)
   s
end

#=
     find the smallest N such that stringTrimmed(lowerbound(x), N) == stringTrimmed(upperbound(x), N)
=#

function smartarbstring{P}(x::ArbFloat{P})
     digits = floor(Int, precision(x)*0.3010299956639811952137)

     if isexact(x)
        return String(x, digits, UInt(2))
     end

     lb, ub = bounds(x)
     lbs = String(lb, digits, UInt(2))
     ubs = String(ub, digits, UInt(2))
     if lbs[end]==ubs[end] && lbs==ubs
         return ubs
     end
     for i in (digits-2):-2:4
         lbs = String(lb, i, UInt(2))
         ubs = String(ub, i, UInt(2))
         if lbs[end]==ubs[end] && lbs==ubs # tests rounding to every other digit position
            us = String(ub, i+1, UInt(2))
            ls = String(lb, i+1, UInt(2))
            if us[end] == ls[end] && us==ls # tests rounding to every digit position
               ubs = lbs = us
            end
            break
         end
     end
     if lbs != ubs
        ubs = String(x, 3, UInt(2))
     end
     ubs
end

function smartvalue{P}(x::ArbFloat{P})
    s = smartarbstring(x)
    ArbFloat{P}(s)
end

function smartstring{P}(x::ArbFloat{P})
    s = smartarbstring(x)
    a = ArbFloat{P}(s)
    postfix =
        if (upperbound(x) < a)
            "-"
        elseif (lowerbound(x) > a)
            "+"
        else
            "~"
        end
    string(s,postfix)
end

function stringAll{P}(x::ArbFloat{P})
    string(midpoint(x)," ± ", string(radius(x),17))
end

function stringCompact{P}(x::ArbFloat{P})
    string(x,8)
end

function stringAllCompact{P}(x::ArbFloat{P})
    string(string(midpoint(x),8)," ± ", string(radius(x),5))
end


text:
using Mass2
using Base.Test

include("runningvectors.jl")
include("mockpulses.jl")
include("calibration.jl")
include("histogram.jl")
include("corrections.jl")
include("ljh.jl")
include("tes_optimal_filtering.jl")
include("steps_macros.jl")
#include("steps_with_mockpulses.jl")
include("good_mnka_mystery.jl")
#include("ljh_calibronium.jl")
include("summarize.jl")


text:
#Horizontal line
function get_line(x0,x1 , E)

    x =  collect(x0:0.1:x1)
    y = similar(x)
    y .= E

    return x,y
end

#Exponential decay
function get_line(x0 , x1 , a , E  )

    x = collect(x0:0.1:x1)
    y = E.*exp.(-a.*(x.-x0))

    return x,y

end

function get_simple_theoretical_model_data(p::Param , tmax = p.tmax)

    #Storage Vectors
    t_vec = Float64[]
    y_vec = Float64[]

    #If no dissipation
    if(!p.dissipation)

        #Initial configurations
        E = 1.0
        t = 0

        t_vec , y_vec = get_line(t , tmax  , E)

    elseif(p.R_orbit == 0.0)

        #Initial configurations
        E = 1.0
        t = 0
        absorb = false

        #Initial pulse configuration
        Ro = 0.9*p.xmax
        Ri = p.R1
        Rorbit = p.R_orbit

        #First iteration
        x , y = get_line(t , t + (p.r0-Ri)  , E)

        append!(t_vec, x)
        append!(y_vec, y)

        #push!(t_vec , x)
        #push!(y_vec , y)

        absorb = true
        t+= (p.r0-Ri)

        while(t < tmax)

            #absorption
            if(absorb)
                x , y = get_line(t , t+2*Ri , p.alpha1 ,E)

                append!(t_vec, x)
                append!(y_vec, y)

                                #push!(t_vec , x)
                #push!(y_vec , y)

                E = E*exp(- p.alpha1*(2*Ri))
                absorb = false
                t+= 2*Ri
            else
                x , y = get_line(t , t+ 2*(Ro-Ri) , E)

                append!(t_vec, x)
                append!(y_vec, y)

                absorb = true
                t+= 2*(Ro-Ri)
            end

        end

    else

        #Initial configuration
        E = 1.0
        t = 0

        Ro = 0.9*p.xmax
        Ri = p.R1
        Rorbit = p.R_orbit

        #Geometrical factor and phase of interaction
        phase = 1
        gfactor =  (2/pi)*(Ri/(Rorbit+Ri))

        alpha = gfactor*p.alpha1

        #First iteration
        x , y = get_line(t , t + (p.r0-Rorbit)  , E)

        push!(t_vec , x)
        push!(y_vec , y)

        phase = 2
        t+= (p.r0-Rorbit)

        while(t < tmax)

            #absorption
            if(phase == 1)

                x , y = get_line(t , t+(Ro-Rorbit-Ri) , E)

                push!(t_vec , x)
                push!(y_vec , y)

                phase = 2
                t+= (Ro-Rorbit-Ri)

            elseif(phase == 2)

                x , y = get_line(t , t+ 2*Ri , alpha , E)

                push!(t_vec , x)
                push!(y_vec , y)

                E = E*exp(- alpha*(2*Ri))
                phase = 3
                t+= 2*Ri

            elseif(phase == 3)

                x , y = get_line(t , t+2*(Rorbit-Ri) , E)

                push!(t_vec , x)
                push!(y_vec , y)

                phase = 4
                t+= 2*(Rorbit-Ri)

            elseif(phase == 4)

                x , y = get_line(t , t+ 2*Ri , alpha , E)

                push!(t_vec , x)
                push!(y_vec , y)

                E = E*exp(- alpha*(2*Ri))
                phase = 5
                t+= 2*Ri

            elseif(phase == 5)

                x , y = get_line(t , t + (Ro-Rorbit+ Ri) ,  E)

                push!(t_vec , x)
                push!(y_vec , y)

                phase = 1
                t+= (Ro-Rorbit)

            end
        end
    end

    return t_vec , y_vec

end


text:
using Test, aoc2019

@testset "juliaify_address" begin
    @test juliaify_address(1) == 2
end

@testset "process_intcode" begin
    (intcode = [1,9,10,3,2,3,11,0,99,30,40,50]; process_intcode!(intcode))
    @test intcode == [3500,9,10,70,2,3,11,0,99,30,40,50]

    (intcode = [1,0,0,0,99]; process_intcode!(intcode))
    @test intcode == [2,0,0,0,99]

    (intcode = [2,3,0,3,99]; process_intcode!(intcode))
    @test intcode == [2,3,0,6,99]

    (intcode = [2,4,4,5,99,0]; process_intcode!(intcode))
    @test intcode == [2,4,4,5,99,9801]

    (intcode = [1,1,1,4,99,5,6,0,99]; process_intcode!(intcode))
    @test intcode == [30,1,1,4,2,5,6,0,99]

    @test_throws BoundsError process_intcode!([1,0,0,0])
    @test_throws BoundsError process_intcode!([1,0,0,0,1])
    @test_throws BoundsError process_intcode!([1,0,0,10,99])
    @test_throws ArgumentError process_intcode!([3,0,0,0,1,2,2,3])
end

@testset "elf_intcode_input" begin
    (intcode = [1,9,10,3,2]; index=0; fix=2; elf_intcode_input!(intcode, index, fix))
    @test intcode == [2,9,10,3,2]
end

@testset "determine_elf_intcode_inputs" begin
    intcode = aoc2019.d02_inputread(aoc2019.d02_inputpath)
    @test determine_elf_intcode_inputs(intcode, 4690667) == (12, 2)
end

@testset "d02_answers" begin
    @test aoc2019.d02_inputpath isa String
    input = aoc2019.d02_inputread(aoc2019.d02_inputpath)
    @test length(input) == 173
    @test aoc2019.d02_answer1(input) == 4690667
    @test aoc2019.d02_answer2(input) == 6255
end


text:
# constructor
function SemML(
        observed::O where {O <: SemObs},
        imply::I where {I <: ImplyDefinition},
        objective,
        grad)
    n_obs = Int64(observed.n_obs)
    n_man = Int64(observed.n_man)
    n_patterns = imply.n_patterns

    inverses = [zeros(n_man, n_man) for i = 1:n_patterns]
    choleskys = Array{Cholesky{Float64,Array{Float64,2}},1}(undef, n_patterns)
    logdets = zeros(n_patterns)

    meandiff = [zeros(n_man) for i = 1:n_obs]

    return SemML(
    inverses,
    choleskys,
    nothing,
    logdets,
    meandiff,
    copy(objective),
    copy(grad)
    )
end

# loss
function (semml::SemML)(par, model::Sem{O, I, L, D}) where
    {O <: SemObs, L <: Loss, I <: ImplyDefinition, D <: SemFiniteDiff}

    if !prepare_ML!(semml, model) return Inf end

    F = zero(eltype(par))
    F = F_Definition(F, semml, model)

    return F
end

function batch_cholesky!(semml::SemML, model)
    for i = 1:size(semml.choleskys, 1)
        semml.choleskys[i] =
            cholesky!(Hermitian(model.imply.imp_cov[i]); check = false)
        if !isposdef(semml.choleskys[i]) return false end
    end
    return true
end

function prepare_ML!(semml, model)
    if !batch_cholesky!(semml, model) return false end
    semml.logdets .= logdet.(semml.choleskys)
    batch_inv!(semml, model)
    return true
end

function F_Definition(F, semml, model)
    for i = 1:Int64(model.imply.n_patterns)
        for j in model.imply.rows[i]
            let (imp_mean, meandiff, inverse, data, logdet) =
                (model.imply.imp_mean[i],
                semml.meandiff[i],
                semml.inverses[i],
                model.observed.data_rowwise[j],
                semml.logdets[i])

                F += F_one_person(imp_mean, meandiff, inverse, data, logdet)

            end
        end
    end
    return F
end

function findrow(r, rows)
    for i in 1:length(rows)
        if r ∈ rows[i]
            return i
        end
    end
end


text:
basename_type(x) = basename_type(typeof(x))
basename_type(::Type{P}) where {D,B,P<:NeuroPath{D,B}} = B

dirname_type(x) = dirname_type(typeof(x))
dirname_type(::Type{P}) where {D,B,P<:NeuroPath{D,B}} = D

similar_type(::Type{<:StudyPath}, ::Type{D}, ::Type{B}) where {D,B} = StudyPath{D,B}
similar_type(::Type{<:DerivativePath}, ::Type{D}, ::Type{B}) where {D,B} = DerivativePath{D,B}
similar_type(::Type{<:SessionPath}, ::Type{D}, ::Type{B}) where {D,B} = SessionPath{D,B}
similar_type(::Type{<:DataPath}, ::Type{D}, ::Type{B}) where {D,B} = DataPath{D,B}
similar_type(::Type{<:SubjectPath}, ::Type{D}, ::Type{B}) where {D,B} = SubjectPath{D,B}
similar_type(::Type{<:ModalityPath}, ::Type{D}, ::Type{B}) where {D,B} = ModalityPath{D,B}
similar_type(::Type{<:FilePath}, ::Type{D}, ::Type{B}) where {D,B} = FilePath{D,B}

# TODO """ ModalityIndicator """
Base.@kwdef struct ModalityIndicator <: AbstractVector{Modality}
    anat::Bool=false
    beh::Bool=false
    dwi::Bool=false
    fmap::Bool=false
    func::Bool=false
    perf::Bool=false
    pet::Bool=false
    eeg::Bool=false
    ieeg::Bool=false
    meg::Bool=false
end
@propagate_inbounds function Base.getindex(x::ModalityIndicator, i::Int)::Modality
    out = _get_indicator(x, i)
    @boundscheck out === nothing && throw(BoundsError(x, i))
    return getfield(out, 1)
end
function _get_indicator(x::ModalityIndicator, i::Int)
    if i === 1
        if getfield(x, i)
            return ModalityPath(nothing, "anat"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 2
        if getfield(x, i)
            return ModalityPath(nothing, "beh"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 3
        if getfield(x, i)
            return ModalityPath(nothing, "dwi"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 4
        if getfield(x, i)
            return ModalityPath(nothing, "fmap"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 5
        if getfield(x, i)
            return ModalityPath(nothing, "func"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 6
        if getfield(x, i)
            return ModalityPath(nothing, "perf"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 7
        if getfield(x, i)
            return ModalityPath(nothing, "pet"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 8
        if getfield(x, i)
            return ModalityPath(nothing, "eeg"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 9
        if getfield(x, i)
            return ModalityPath(nothing, "ieeg"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 10
        if getfield(x, i)
            return ModalityPath(nothing, "meg"), i
        else
            return _get_indicator(x, i+1)
        end
    else
        return nothing
    end
end

# TODO """ DerivativeIndicator """
Base.@kwdef struct DerivativeIndicator <: AbstractVector{Data}
    rawdata::Bool=false
    sourcedata::Bool=false
    stimuli::Bool=false
    derivatives::Bool=false
    phenotype::Bool=false
    code::Bool=false
end
@propagate_inbounds function Base.getindex(x::DerivativeIndicator, i::Int)::Data
    out = _get_indicator(x, i)
    @boundscheck out === nothing && throw(BoundsError(x, i))
    return getfield(out, 1)
end
function _get_indicator(x::DerivativeIndicator, i::Int)
    if i === 1
        if getfield(x, i)
            return DataPath(nothing, "rawdata"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 2
        if getfield(x, i)
            return DataPath(nothing, "sourcedata"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 3
        if getfield(x, i)
            return DataPath(nothing, "stimuli"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 4
        if getfield(x, i)
            return DataPath(nothing, "derivatives"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 5
        if getfield(x, i)
            return DataPath(nothing, "phenotype"), i
        else
            return _get_indicator(x, i+1)
        end
    elseif i === 6
        if getfield(x, i)
            return DataPath(nothing, "code"), i
        else
            return _get_indicator(x, i+1)
        end
    else
        return nothing
    end
end

function Base.iterate(x::Union{DerivativeIndicator,ModalityIndicator}, state=0)
    _get_indicator(x, state + 1)
end

# TODO """ PathVector """
struct PathVector{T,D,B} <: AbstractVector{T}
    dirname::D
    basenames::B

    function PathVector(dirname::D, basenames::AbstractVector{T}) where {D,T}
        new{similar_type(T, D, basename_type(T)),D,typeof(B)}(dirname, basenames)
    end
end

@inline function Base.iterate(x::PathVector)
    itr = iterate(getfield(x, :basenames))
    if itr === nothing
        return nothing
    else
        return (joinpath(dirname(x), getfield(itr, 1)), getfield(itr, 2))
    end
end
@inline function Base.iterate(x::PathVector, state)
    itr = iterate(getfield(x, :basenames), state)
    if itr === nothing
        return nothing
    else
        return (joinpath(dirname(x), getfield(itr, 1)), getfield(itr, 2))
    end
end

@propagate_inbounds function Base.getindex(x::PathVector, i::AbstractVector)
    @boundscheck checkbounds(getfield(x, :basenames), i)
    joinpath(dirname(x), @inbounds(getfield(x, :basenames)[i]))
end

@propagate_inbounds function Base.getindex(x::PathVector, i::Int)
    @boundscheck checkbounds(getfield(x, :basenames), i)
    joinpath(dirname(x), @inbounds(getfield(x, :basenames)[i]))
end


text:
module NeutralEvolution
include("neutral.jl")
include("propsel_funct.jl")
end
using Main.NeutralEvolution


text:
module Logical

using StatsFuns: logsumexp

using ..Data
using ..Utils

export

# Circuits
Var, Lit, var2lit, lit2var, ΔNode, Δ,
LiteralGate, ConstantGate, ⋁Gate, ⋀Gate,
GateType, InnerGate, LeafGate,
negative, positive, variable, literal,
true_like, false_like, DecoratorΔ, DecoratorΔNode, origin, grand_origin, is_true, is_false, constant,
prime, sub, conjoin_like, disjoin_like, literal_like,

# CircuitTraversal
⋁_nodes, ⋀_nodes,
or_nodes, and_nodes,

# Queries
num_variables, is_decomposable, variable_scopes, variable_scope,
model_count, sat_prob, prob_equiv_signature, has_unique_literal_nodes, has_unique_constant_nodes,
is_smooth,

# Transformations
smooth, forget, propagate_constants, condition, split, clone, replace_node,

# LogicalCircuits
LogicalΔNode, UnstLogicalΔNode, LogicalLeafNode, LogicalΔ, UnstLogicalΔ, node_type,
LiteralNode, ConstantNode, ⋁Node, ⋀Node, TrueNode, FalseNode, fully_factorized_circuit, copy,
copy_node, normalize,

# UpFlowCircuits
UpFlowΔNode, UpFlowΔ, UpFlowLeafNode, UpFlowInnerNode,
UpFlowLiteral, UpFlowConstant, UpFlow⋀Compact, UpFlow⋀Cached, UpFlow⋁Compact, UpFlow⋁Cached, UpFlow⋀, UpFlow⋁,
flow_opts★, pass_up, pass_up_node, pr_factors, pr, resize_flows, flow_length,

# DownFlowCircuits
DownFlowΔNode, DownFlowΔ, DownFlowΔ, DownFlowLeaf, DownFlowInnerNode,
DownFlow⋀Compact, DownFlow⋀Cached, DownFlow⋁Compact, DownFlow⋁Cached, DownFlow⋀, DownFlow⋁,
HasDownFlow, pass_down, pass_up_down, downflow, reset_downflow_in_progress, downflow_sinks,
FlowΔNode, FlowΔ, FlowΔ,

# AggregateFlowCircuits
AggregateFlowΔ, AggregateFlowΔ, AggregateFlowΔNode, AggregateFlow⋁,
reset_aggregate_flows, accumulate_aggr_flows_batch, accumulate_aggr_flows_cached, opts_accumulate_flows, collect_aggr_flows,
accumulate_aggr_flows,

# Vtree
VtreeNode, Vtree, balanced_vtree, top_down_vtree, bottom_up_vtree,

# StructuredLogicalCircuits
StructLogicalΔNode, StructLogicalLeafNode, StructLogicalΔ,
StructLiteralNode, StructConstantNode, Struct⋁Node, Struct⋀Node, vtree,

# PlainVtree
PlainVtreeNode, PlainVtreeLeafNode, PlainVtreeInnerNode, isleaf, variables, num_variables, PlainVtree,
path_length, random_vtree,

# Sdd
SddMgrNode, SddMgr,
SddNode, SddLeafNode, SddCircuit, Sdd,
SddLiteralNode, SddConstantNode, Sdd⋁Node, Sdd⋀Node,

# SddMgr
sdd_size, sdd_num_nodes,
compile_cnf, compile_clause, validate,

# TrimSddMgr
TrimMgrNode, TrimSddMgr, XYPartition, Element, TrimSdd, TrimNode,
compress, unique⋁, canonicalize,
compile, conjoin, disjoin, negate, descends_from, descends_left_from, descends_right_from,

# temp
simple_test

include("Circuits.jl")
include("CircuitTraversal.jl")

include("LogicalCircuits.jl")
include("UpFlowCircuits.jl")
include("DownFlowCircuits.jl")
include("AggregateFlowCircuits.jl")

include("Queries.jl")
include("Transformations.jl")

include("vtrees/Vtree.jl")
include("StructuredLogicalCircuits.jl")

include("vtrees/PlainVtree.jl")
include("Sdd.jl")
include("vtrees/SddMgr.jl")
include("vtrees/TrimSddMgr.jl")

end


text:
function generateIntegers(n, min, max,replacement)
    url = "https://api.random.org/json-rpc/2/invoke"
    apikey = "f72f3135-206b-4f61-adc0-143012355223"
    headers = (("content-type", "application/json"),
                       ("Accept", "application/json"))

           params = Dict("apiKey"=>apikey, "n"=>n, "min"=>min, "max"=>max, "replacement"=>replacement)
           req = Dict("jsonrpc" => "2.0", "method" =>"generateIntegers", "params" => params, "id"=>0)
           response = HTTP.post(url, headers, JSON.json(req))
           response_json = String(response.body)
           response_parsed = JSON.parse(response_json)
           random_op = response_parsed["result"]["random"]["data"]
    return(random_op)
end

function generateIntegerSequences(n,length, min, max,replacement)
    url = "https://api.random.org/json-rpc/2/invoke"
    apikey = "f72f3135-206b-4f61-adc0-143012355223"
    headers = (("content-type", "application/json"),
                       ("Accept", "application/json"))

           params = Dict("apiKey"=>apikey, "n"=>n, "length"=>length, "min"=>min, "max"=>max, "replacement"=>replacement)
           req = Dict("jsonrpc" => "2.0", "method" =>"generateIntegerSequences", "params" => params, "id"=>0)
           response = HTTP.post(url, headers, JSON.json(req))
           response_json = String(response.body)
           response_parsed = JSON.parse(response_json)
           random_op = response_parsed["result"]["random"]["data"]
    return(random_op)
end

function generateDecimalFractions(n, decimalPlaces, replacement)
    url = "https://api.random.org/json-rpc/2/invoke"
    apikey = "f72f3135-206b-4f61-adc0-143012355223"
    headers = (("content-type", "application/json"),
                       ("Accept", "application/json"))

           params = Dict("apiKey"=>apikey, "n"=>n, "decimalPlaces"=>decimalPlaces, "replacement"=>replacement)
           req = Dict("jsonrpc" => "2.0", "method" => "generateDecimalFractions", "params" => params, "id"=>0)
           response = HTTP.post(url, headers, JSON.json(req))
           response_json = String(response.body)
           response_parsed = JSON.parse(response_json)
           random_op = response_parsed["result"]["random"]["data"]
    return(random_op)
end

function generateGaussians(n, mean, standard_dev, significant_digits)
    url = "https://api.random.org/json-rpc/2/invoke"
    apikey = "f72f3135-206b-4f61-adc0-143012355223"
    headers = (("content-type", "application/json"),
                       ("Accept", "application/json"))

           params = Dict("apiKey"=>apikey, "n"=>n, "mean"=>mean, "standardDeviation"=>standard_dev, "significantDigits"=>significant_digits)
           req = Dict("jsonrpc" => "2.0", "method" =>"generateGaussians", "params" => params, "id"=>0)
           response = HTTP.post(url, headers, JSON.json(req))
           response_json = String(response.body)
           response_parsed = JSON.parse(response_json)
           random_op = response_parsed["result"]["random"]["data"]
    return(random_op)
end

function generateStrings(n, length, characters, replacement)
    url = "https://api.random.org/json-rpc/2/invoke"
    apikey = "f72f3135-206b-4f61-adc0-143012355223"
    headers = (("content-type", "application/json"),
                       ("Accept", "application/json"))

           params = Dict("apiKey"=>apikey, "n"=>n, "length"=>length, "characters"=>characters, "replacement"=>replacement)
           req = Dict("jsonrpc" => "2.0", "method" =>"generateStrings", "params" => params, "id"=>0)
           response = HTTP.post(url, headers, JSON.json(req))
           response_json = String(response.body)
           response_parsed = JSON.parse(response_json)
           random_op = response_parsed["result"]["random"]["data"]
    return(random_op)
end

function generateUUIDs(n)
    url = "https://api.random.org/json-rpc/2/invoke"
    apikey = "f72f3135-206b-4f61-adc0-143012355223"
    headers = (("content-type", "application/json"),
                       ("Accept", "application/json"))

           params = Dict("apiKey"=>apikey, "n"=>n)
           req = Dict("jsonrpc" => "2.0", "method" =>"generateUUIDs", "params" => params, "id"=>0)
           response = HTTP.post(url, headers, JSON.json(req))
           response_json = String(response.body)
           response_parsed = JSON.parse(response_json)
           random_op = response_parsed["result"]["random"]["data"]
    return(random_op)
end

function generateBlobs(n, size_in_bits, format)
    url = "https://api.random.org/json-rpc/2/invoke"
    apikey = "f72f3135-206b-4f61-adc0-143012355223"
    headers = (("content-type", "application/json"),
                       ("Accept", "application/json"))

           params = Dict("apiKey"=>apikey, "n"=>n, "size"=>size_in_bits, "format"=>format)
           req = Dict("jsonrpc" => "2.0", "method" =>"generateBlobs", "params" => params, "id"=>0)
           response = HTTP.post(url, headers, JSON.json(req))
           response_json = String(response.body)
           response_parsed = JSON.parse(response_json)
           random_op = response_parsed["result"]["random"]["data"]
    return(random_op)
end