text:
# functions to read/write networks topologies
# originally in functions.jl
# Claudia March 2015

# ------------------------------- Read topology ------------------------------------------

function peekchar(s::IOStream)
    Base.peekchar(s)
end

function peekchar(s::IOBuffer)
    mark(s)
    c = read(s,Char)
    reset(s)
    return c
end

# aux function to detect characters which should be ignored by readTopology
function iswhitesymbol(c::Char)
    whitesymbols = [' ', '\n', '\r', '\t']
    return c in whitesymbols
end

# aux function to peek the next non-white-symbol char in s, does not modify s
# input: s IOStream/IOBuffer
function peekskip(s::IO)
    c = peekchar(s)
    mark(s)
    while iswhitesymbol(c)
        c = read(s, Char)
    end
    reset(s)
	return c
end

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

# aux function to advance stream in readSubtree
# input: s IOStream/IOBuffer
function advance!(s::IO, c::Char, numLeft::Array{Int,1})
    c = peekskip(s)
    if(Base.eof(s))
        error("Tree ends prematurely while reading subtree after left parenthesis $(numLeft[1]-1).")
    end
    return readskip!(s)
end

# aux function to read all digits of taxon name
# it allows names with letters and numbers
# it also reads # as part of the name and returns pound=true
# it returns the node name as string as well to check if it exists already (as hybrid)
# warning: treats digit taxon numbers as strings to avoid repeated node numbers
function readNum(s::IO, c::Char, net::HybridNetwork, numLeft::Array{Int,1})
    pound = 0;
    if(isalnum(c) || isValidSymbol(c) || c == '#')
        pound += (c == '#') ? 1 : 0
        name = readskip!(s)
        c = peekskip(s)
        while(isalnum(c) || isValidSymbol(c) || c == '#')
            d = readskip!(s)
            name = string(name,d)
            if(d == '#')
                pound += 1;
               c = peekskip(s);
               if(isalnum(c))
                   if(c != 'L' && c != 'H' && c != 'R')
                       warn("Expected H, R or LGT after # but received $(c) in left parenthesis $(numLeft[1]-1).")
                   end
               else
                   a = readstring(s);
                   error("Expected name after # but received $(c) in left parenthesis $(numLeft[1]-1). Remaining is $(a).")
               end
            end
            c = peekskip(s);
        end
        if(pound == 0)
            return size(net.names,1)+1, name, false
        elseif(pound == 1)
            return size(net.names,1)+1, name, true
        else
            a = readstring(s);
            error("strange node name with $(pound) # signs. remaining is $(a).")
        end
    else
        a = readstring(s);
        error("Expected int digit, alphanum or # but received $(c). remaining is $(a).");
    end
end

# aux function to read floats like length
function readFloat(s::IO, c::Char)
    if(isdigit(c) || in(c, ['.','e','-']))
        num = string(readskip!(s));
        c = peekskip(s);
        while(isdigit(c) || in(c, ['.','e','-']))
            d = readskip!(s);
            num = string(num,d);
            c = peekskip(s);
        end
        f = 0.0
        try
            f = float(num)
        catch
            error("problem with number read $(num), not a float number")
        end
        return f
    else
        a = readstring(s);
        error("Expected float digit after : but found $(c). remaining is $(a).");
    end
end

# aux function to identify if a symbol in taxon name is valid
# symbol cannot be: () [] : ; ' , . space \t \r \n
# according to richnewick.pdf
function isValidSymbol(c::Char)
    return !isspace(c) && c != '(' && c != ')' && c != '[' && c != ']' && c != ':' && c != ';' && c != ',' #&& c != '.'
end

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

Helper function for readSubtree!
Called once a `(` has been read in a tree topology and reads until the corresponding `)` has been found.
This function performs the recursive step for readSubtree!
Advances `s` past the subtree, adds discovered nodes and edges to `net`, `hybrids`, and `index`.
"""
@inline function parseRemainingSubtree!(s::IO, numLeft::Array{Int,1}, net::HybridNetwork, hybrids::Array{String, 1}, index::Array{Int, 1})
    numLeft[1] += 1
    DEBUGC && println(numLeft)
    n = Node(-1*numLeft[1],false);
    DEBUG && println("creating node $(n.number)")
    keepon = true;
    c = readskip!(s)
    while (keepon)
        bl = readSubtree!(s,n,numLeft,net,hybrids,index)
        c = advance!(s,c,numLeft)
        if (c == ')')
            keepon = false
        elseif (c != ',')
            a = readstring(s);
            error("Expected right parenthesis after left parenthesis $(numLeft[1]-1) but read $(c). The remainder of line is $(a).")
        end
    end
    return n
end

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

Helper function for readSubtree!
Handles any type of given hybrid node.
To be called once a `#` has been found in a tree topology.
Creates a hybrid node and its edges, then inserts those into `net`, `hybrids` and `index` accordingly.
"""
@inline function parseHybridNode!(n::Node, parent::Node, name::String, net::HybridNetwork, hybrids::Array{String, 1}, index::Array{Int, 1})
    DEBUG && println("found pound in $(name)")
    n.hybrid = true;
    DEBUGC && println("got hybrid $(name)")
    DEBUGC && println("hybrids list has length $(length(hybrids))")
    ind = findfirst(hybrids, name) # index of 'name' in the list 'hybrid'. 0 if not found
    if ind>0 # the hybrid name was seen before
        DEBUG && println("$(name) was found in hybrids list")
        other = net.node[index[ind]]; # other node with same hybrid name
        DEBUG && println("other is $(other.number)")
        DEBUGC && println("other is leaf? $(other.leaf), n is leaf? $(n.leaf)")
        if !n.leaf && !other.leaf
            error("both hybrid nodes are internal nodes: successors of the hybrid node must only be included in the node list of a single occurrence of the hybrid node.")
        elseif n.leaf
            DEBUG && println("n is leaf")
            e = Edge(net.numEdges+1);
            DEBUG && println("creating hybrid edge $(e.number) attached to other $(other.number) and parent $(parent.number)")
            e.hybrid = true
            e.isMajor = false; # isMajor = true by default constructor
            pushEdge!(net,e);
            setNode!(e,[other,parent]); # isChild1 = true by default constructor
            setEdge!(other,e);
            setEdge!(parent,e);
            DEBUG && println("e $(e.number )istIdentifiable? $(e.istIdentifiable)")
        else # !n.leaf
            DEBUG && println("n is not leaf, other is leaf")
            size(other.edge,1) == 1 || # other should be a leaf
               error("strange: node $(other.number) is a leaf hybrid node. should have only 1 edge but has $(size(other.edge,1))")
            DEBUGC && println("other is $(other.number), n is $(n.number), edge of other is $(other.edge[1].number)")
            otheredge = other.edge[1];
            otherparent = getOtherNode(otheredge,other);
            DEBUG && println("otheredge is $(otheredge.number)")
            DEBUG && println("parent of other is $(otherparent.number)")
            removeNode!(other,otheredge);
            deleteNode!(net,other);
            setNode!(otheredge,n);
            setEdge!(n,otheredge);
            ## otheredge.istIdentifiable = true ## setNode should catch this, but when fixed, causes a lot of problems
            DEBUG && println("setting otheredge to n $(n.number)")
            e = Edge(net.numEdges+1) # isMajor = true by default
            DEBUG && println("creating hybrid edge $(e.number) between n $(n.number) and parent $(parent.number)")
            e.hybrid = true
            setNode!(e,[n,parent]);
            setEdge!(n,e);
            setEdge!(parent,e);
            pushNode!(net,n);
            pushEdge!(net,e);
            n.number = other.number;
            n.name = other.name;
            DEBUG && println("edge $(e.number) istIdentifiable? $(e.istIdentifiable)")
            DEBUG && println("otheredge $(otheredge.number) istIdentifiable? $(otheredge.istIdentifiable)")
        end
    else # ind=0: hybrid name not seen before
        DEBUG && println("$(name) not found in hybrids list")
        DEBUG && println("$(name) is leaf? $(n.leaf)")
        n.hybrid = true;
        nam = string(name)
        push!(net.names, nam);
        n.name = nam;
        DEBUGC && println("put $(nam) in hybrids name list")
        push!(hybrids, nam);
        pushNode!(net,n);
        push!(index,size(net.node,1));
        e = Edge(net.numEdges+1);
        DEBUG && println("creating hybrid edge $(e.number)")
        e.hybrid = true
        n.leaf ? e.isMajor = false : e.isMajor = true
        pushEdge!(net,e);
        setNode!(e,[n,parent]);
        setEdge!(n,e);
        setEdge!(parent,e);
        DEBUG && println("edge $(e.number) istIdentifiable? $(e.istIdentifiable)")
    end
    e.containRoot = !e.hybrid # not good: but necessay for SNaQ functions
    return e
end

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

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

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

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

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

#TODO: Check the windows version

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

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

    edgeLenPresent = (edgeLen != -1.0);
    gammaPresent   = (gamma != -1.0 && gamma != nothing);

    e.length = edgeLen
    if gammaPresent
        if(!e.hybrid)
            warn("gamma read for current edge $(e.number) but it is not hybrid, so gamma=$(gamma) ignored")
        else
            parentEdges = Edge[] # The edges which have n as a child
            for e in n.edge
                if e.hybrid
                    nIsChild = (e.isChild1 && e.node[1] == n || !e.isChild1 && e.node[2] == n)
                    if nIsChild
                        push!(parentEdges, e)
                    end
                end
            end
            numParents = size(parentEdges)[1]
            if (numParents == 0) || (numParents == 1)
                e.gamma = gamma
                # isMajor will be set when the other edge is found
            elseif (numParents == 2) #other edge has been read, may have gamma set
                local otheredge::Edge
                parentEdges[1] == e ? otheredge = parentaledges[2] :
                                        otheredge = parentEdges[1]
                # Note: only need to correct isMajor if the gamma values are not equal
                if (!approxEq(gamma + otheredge.gamma, 1)) # gammas do not sum to 1
                    # some gamma < 0
                    if (gamma <= 0 && otheredge.gamma > 0)
                        gamma = 1 - otheredge.gamma
                    elseif (otheredge.gamma <= 0 && gamma > 0)
                        otheredge.gamma = 1 - gamma
                    end
                    # rescale so that they sum to 1
                    e.gamma, otheredge.gamma = ((gamma)/(gamma + otheredge.gamma),
                                                (otheredge.gamma)/(gamma + otheredge.gamma))
                else # both gamma values sum to 1
                    e.gamma = gamma
                end
                if (!approxEq(e.gamma, 0.5))
                        if e.gamma > 0.5
                            e.isMajor = true
                            otheredge.isMajor = false
                        else
                            e.isMajor = false
                            otheredge.isMajor = true
                        end
                    end
            else
                error("hybrid edge gamma value parsed but hybrid node has $(size(parentEdges)[1]) parent edges! (should be between 0 and 2)")
            end
        end
    end
end

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

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

Reads additional info formatted as: `:length:bootstrap:gamma`.
Allows for name of internal nodes without # after closing parenthesis: (1,2)A.
Warning if hybrid edge without γ, or if γ (ignored) without hybrid edge
"""
# modified from original Cecile c++ code to allow polytomies

function readSubtree!(s::IO, parent::Node, numLeft::Array{Int,1}, net::HybridNetwork, hybrids::Array{String,1}, index::Array{Int,1})
    c = peekskip(s)
    e = nothing;
    hasname = false; # to know if the current node has name
    pound = false;
    if(c == '(')
        # read the rest of the subtree (perform the recursive step!)
        n = parseRemainingSubtree!(s, numLeft, net, hybrids, index)
        c = peekskip(s);
        if(isalnum(c) || isValidSymbol(c) || c == '#') # internal node has name
            hasname = true;
            num, name, pound = readNum(s, c, net, numLeft);
            n.number = num;
            c = peekskip(s);
            #if(!pound)
            #    warn("internal node with name without it being a hybrid node. node name might be meaningless after tree modifications.")
            #end
        end
    elseif(isalnum(c) || isValidSymbol(c) || c == '#')
        hasname = true;
        bl = true;
        num, name, pound = readNum(s, c, net, numLeft)
        n = Node(num, true);
        DEBUG && println("creating node $(n.number)")
    else
        a = readstring(s);
        error("Expected beginning of subtree but read $(c) after left parenthesis $(numLeft[1]-1), remaining is $(a).");
    end
    if(pound) # found pound sign in name
        e = parseHybridNode!(n, parent, name, net, hybrids, index)
    else
        if(hasname)
            push!(net.names,string(name));
            n.name = string(name)
        end
        e = parseTreeNode!(n, parent, net)
    end
    c = peekskip(s);
    e.length = -1.0
    e.gamma = e.hybrid? -1.0 : 1.0
    if(c == ':')
        parseEdgeData!(s, e, n, numLeft)
    end
    return true
end

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

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

Read tree or network topology from parenthetical format (extended Newick).

Input: text file or parenthetical format directly.
The file name may not start with a left parenthesis, otherwise the file
name itself would be interpreted as the parenthetical description.
"""
readTopology(input::AbstractString) = readTopology(input,true)

function readTopology(s::IO,verbose::Bool)
    net = HybridNetwork()
    line = readuntil(s,";");
    if(line[end] != ';')
        error("file does not end in ;")
    end
    seekstart(s)
    c = peekskip(s)
    numLeft = [1]; # made Array to make it mutable; start at 1 to avoid node -1 which breaks undirectedOtherNetworks
    hybrids = String[];
    index = Int[];
    if(c == '(')
        numLeft[1] += 1;
        #println(numLeft)
        n = Node(-1*numLeft[1],false);
        c = readskip!(s)
        b = false;
        while(c != ';')
            b |= readSubtree!(s,n,numLeft,net,hybrids,index)
            c = readskip!(s);
            if(eof(s))
                error("Tree ended while reading in subtree beginning with left parenthesis number $(numLeft[1]-1).")
            elseif(c == ',')
                continue;
            elseif(c == ')')
                c = peekskip(s);
                if(c == ':')
                    while(c != ';')
                        c = readskip!(s)
                    end
                end
            end
        end
        DEBUG && println("after readsubtree:")
        DEBUG && printEdges(net)
        if(size(n.edge,1) == 1) # root has only one child
            edge = n.edge[1]; # assume it has only one edge
            child = getOtherNode(edge,n);
            removeEdge!(child,edge);
            net.root = getIndex(child,net);
            deleteEdge!(net,edge);
        else
            pushNode!(net,n);
            net.root = getIndex(n,net);
        end
    else
		a = readstring(s)
        error("Expected beginning of tree with ( but received $(c) instead, rest is $(a)")
    end
    storeHybrids!(net)
    checkNumHybEdges!(net)
       ## if(verbose)
       ##     any([(e.length == -1.0 && e.istIdentifiable) for e in net.edge]) && println("edges lengths missing, so assigned default value of -1.0") #fixit: not best approach, better to add a flag inside readSubTree, careful bool not modified inside function, need bool array
       ## end
    net.isRooted = true
    return net
end

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

"""
    `checkNumHybEdges!(net)`

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

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

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

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

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

# function to clean topology after readTopology
# looks for:
# TREE:
# - all tree edges must have gamma=1. fixit: cannot point out which doesn't,
#   only shows error.
# - internal nodes with only 2 edges and solves this case.
# - polytomies and choose one resolution at random, issuing a warning
# NETWORK:
# - number of hybrid edges per hybrid node:
#   if 0,1: error (with warning in old functions)
#   if >2: error of hybrid polytomy
#   if 2: check number of tree edges
# - number of tree edges per hybrid node:
#   if 0: leaf hybrid, add child
#   if >1: expand child
#   if 1: check values of gamma:
# - gammas: need to sum to one and be present.
#   error if they do not sum up to one
#   default values of 0.1,0.9 if not present
# leaveRoot=true: leaves the root even if it has only 2 edges (for plotting), default=false
function cleanAfterRead!(net::HybridNetwork, leaveRoot::Bool)
    mod(sum([!e.hybrid?e.gamma:0 for e in net.edge]),1) == 0 ? nothing : error("tree (not network) read and some tree edge has gamma different than 1")
    nodes = copy(net.node)
    for n in nodes
        if(isNodeNumIn(n,net.node)) # very important to check
            if(size(n.edge,1) == 2)
                if(!n.hybrid)
                    if(!leaveRoot || !isEqual(net.node[net.root],n)) #if n is the root
                        deleteIntNode!(net,n);
                    end
                else
                    hyb = sum([e.hybrid?1:0 for e in n.edge]);
                    if(hyb == 1)
                        deleteIntNode!(net,n);
                    end
                end
            end
            if(!n.hybrid)
                if(size(n.edge,1) > 3)
                    DEBUG && warn("polytomy found in node $(n.number), random resolution chosen")
                    solvePolytomy!(net,n);
                end
                hyb = sum([e.hybrid?1:0 for e in n.edge]);
                if(hyb == 1)
                    n.hasHybEdge == true;
                elseif(hyb > 1)
                    warn("strange tree node $(n.number) with more than one hybrid edge, intersecting cycles maybe")
                end
            else
                hyb = sum([e.hybrid?1:0 for e in n.edge]);
                tre = sum([!e.hybrid?1:0 for e in n.edge]);
                if(hyb > 2)
                    error("hybrid node $(n.number) has more than two hybrid edges attached to it: polytomy that cannot be resolved without intersecting cycles.")
                elseif(hyb == 1)
                    hybnodes = sum([n.hybrid?1:0 for n in net.node]);
                    if(hybnodes == 1)
                        error("only one hybrid node number $(n.number) with name $(net.names[n.number]) found with one hybrid edge attached")
                    else
                        error("current hybrid node $(n.number) with name $(net.names[n.number]) has only one hybrid edge attached. there are other $(hybnodes-1) hybrids out there but this one remained unmatched")
                    end
                elseif(hyb == 0)
                    warn("hybrid node $(n.number) is not connected to any hybrid edges, it was transformed to tree node")
                    n.hybrid = false;
                else # 2 hybrid edges
                    if(tre == 0) #hybrid leaf
                        warn("hybrid node $(n.number) is a leaf, so we add an extra child")
                        addChild!(net,n);
                    elseif(tre > 1)
                        warn("hybrid node $(n.number) has more than one child so we need to expand with another node")
                        expandChild!(net,n);
                    end
                    suma = sum([e.hybrid?e.gamma:0 for e in n.edge]);
                    if(suma == -2)
                        #warn("hybrid edges in read network without gammas")
                        println("hybrid edges for hybrid node $(n.number) have missing gamma's, set default: 0.9,0.1")
                        for e in n.edge
                            if(e.hybrid)
                                if (e.isMajor)
                                    e.gamma = 0.9;
                                    e.isMajor = true
                                else
                                    e.gamma = 0.1;
                                    e.isMajor = false
                                end
                            end
                        end
                    elseif(suma != 1)
                        ed1 = nothing
                        ed2 = nothing
                        for e in n.edge
                            if(e.hybrid)
                                isa(ed1,Void) ? ed1=e : ed2=e
                            end
                        end
                        if(ed1.gamma > 0 && ed2.gamma > 0 && ed1.gamma < 1 && ed2.gamma < 1) #both gammas were set, but contradictory
                            error("hybrid edges for hybrid node $(n.number) have gammas that do not sum up to one: $(ed1.gamma),$(ed2.gamma)")
                        elseif(ed1.gamma != -1.0)
                            warn("only one hybrid edge of hybrid node $(n.number) has gamma value $(ed1.gamma) set, the other edge will be assigned $(1-ed1.gamma).")
                            setGamma!(ed2,1-ed1.gamma, false);
                        else
                            warn("only one hybrid edge of hybrid node $(n.number) has gamma value $(ed2.gamma) set, the other edge will be assigned $(1-ed2.gamma).")
                            setGamma!(ed1,1-ed2.gamma, false);
                        end
                    elseif(suma == 1)
                        for e in n.edge
                            if(e.hybrid)
                                if(approxEq(e.gamma,0.5))
                                    e.isMajor = false
                                    break
                                else
                                    break
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

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

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

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

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

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

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

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

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

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

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

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

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

Write to IO the extended newick format (parenthetical description)
of a network.
If written for dendroscope, inheritance γ's are not written.
If `names` is true, taxon names are written, otherwise taxon numbers
are written instead of taxon names.
If unspecified, branch lengths and γ's are rounded to 3 digits.
"""
# method to start at the root and write the whole tree/network
function writeSubTree!(s::IO, net::HybridNetwork, di::Bool, names::Bool,
                       roundBL::Bool, digits::Integer)
    if net.numNodes == 1
        print(s, (names? net.node[net.root].name : string(net.node[net.root].number)))
    elseif net.numNodes > 1
        print(s,"(")
        degree = length(net.node[net.root].edge)
        for e in net.node[net.root].edge
            writeSubTree!(s,getOtherNode(e,net.node[net.root]),e,di,names,roundBL,digits)
            degree -= 1
            degree == 0 || print(s,",")
        end
        print(s,")")
    end
    print(s,";")
    return nothing
end

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

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

If `parent` is `nothing`, the edge attribute `isChild1` is used
and assumed to be correct to write the subtree rooted at `node`.
This is useful to write a subtree starting at a non-root node.
Example:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# function to read the .out file from snaq (optTopRuns) function
function readOutfile(file::AbstractString)
    try
        s = open(file)
    catch
        error("Could not find or open $(file) file");
    end
    s = open(file)
    line = readline(s)
    DEBUG && println("line read $(line)")
    c = line[1]
    if(c == '(')
       println("Estimated network from file $(file): $(line)")
       net = readTopologyUpdate(line)
       vec = split(line,"-Ploglik = ")
       net.loglik = parse(Float64,vec[2])
    else
        error("output file $(filename).out does not contain a tree in the first line, instead it has $(line); or we had trouble finding ploglik.")
    end
    return net
end

"""
`readSnaqNetwork(output file)`

function to read the estimated network from an .out file generated by the snaq function
"""
readSnaqNetwork(file::AbstractString) = readOutfile(file)

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

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

"""
`readMultiTopology(file)`

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

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

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

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

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

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

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

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

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

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

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

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

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

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

"""
    symmetricTree(n, i=1)

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

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

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

"""
    symmetricNet(n, h, gamma)

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


text:
function network_model(
    inputs::SimulationInputs,
    I_balance::AbstractArray{T},
    voltages::AbstractArray{T},
) where {T <: ACCEPTED_REAL_TYPES}
    # This operation might need improvement, when the total shunts aren't added the
    # function is not allocating
    ybus = get_ybus(inputs) .+ get_total_shunts(inputs)
    bus_count = get_bus_count(inputs)
    rows_vals = SparseArrays.rowvals(ybus)
    ybus_vals = SparseArrays.nonzeros(ybus)
    for n in UnitRange{Int}(1, bus_count)
        # The indexing works this way because of the shape of ybus_rectangular and because of the
        # ordering of
        for i in SparseArrays.nzrange(ybus, n)
            I_balance[n] -= ybus_vals[i] * voltages[rows_vals[i]]
        end
        for i in SparseArrays.nzrange(ybus, n + bus_count)
            I_balance[n + bus_count] -= ybus_vals[i] * voltages[rows_vals[i]]
        end
    end
    return I_balance
end

#=
A = sparse(I,J,V)
rows = rowvals(A)
vals = nonzeros(A)
m, n = size(A)
for j = 1:n
   for i in nzrange(A, j)
      row = rows[i]
      val = vals[i]
      # perform sparse wizardry...
   end
end
=#

ext: jl


text:
using Test
using Glowe

@testset "Model training" begin
    include("train.jl")
end

@testset "Model API" begin
    include("model.jl")
end


text:
@testset "util" begin
    trials = 100

    for m = 1:100, n = 1:100
        A = sprand(Bool, m, n, 0.5)
        @test ChainPartitioners.adjointpattern(A) == permutedims(A)
        A = sprand(UInt, m, n, 0.5)
        A.nzval .= true
        @test ChainPartitioners.adjointpattern(A) == permutedims(A)
    end
end


text:
let
    solver = DPWSolver(n_iterations=n_iter, depth=depth, exploration_constant=ec)
    mdp = LegacyGridWorld()

    policy = solve(solver, mdp)

    state = GridWorldState(1,1)

    a = @inferred action(policy, state)

    clear_tree!(policy)
    @test policy.tree == nothing

    # no action pw
    solver = DPWSolver(n_iterations=n_iter, depth=depth, keep_tree=true, exploration_constant=ec, enable_action_pw=false)
    mdp = LegacyGridWorld()

    policy = solve(solver, mdp)

    state = GridWorldState(1,1)

    a = @inferred action(policy, state)
end


text:
using DataFrames
using StatsModels
using LibGEOS
using GaussianProcesses
using Distributions
using Random
using Optim
using LinearAlgebra

function run_simulation()
    Random.seed!(1) # for replicability
    τ = 0.3 # constant treatment effect
    geordd_df, border_X1, border_X2, border_geo = sim_data(1000, τ)
    fmla = @formula(outcome ~ GP(X1, X2) | region + covarA + covarB + categ)
    geordd = GeoRDD.regions_from_dataframe(fmla, geordd_df)
    # choose a Gaussian process kernel from the GaussianProcesses.jl package
    k_se = SEIso(log(0.3), log(1.0)) # Squared Exponential spatial kernel
    k_m = Const(log(20.0))           # constant kernel for the mean parameters
    βkern = LinIso(log(1.0))         # linear kernel for the other covariates
    logNoise = 1.0
    # create the Gaussian process object from the regional data:
    mgpcv = GeoRDD.MultiGPCovars(geordd, k_se+k_m, βkern, logNoise)
    # optimize the hyperparameters
    opt_output=optimize!(
            mgpcv,
            noise=true,
            kern=true,
            domean=false,
            beta=true,
            options=Optim.Options(
                show_trace=false,
                iterations=1000,
                x_tol=1e-8,
                f_tol=1e-10),
        )
    # extract posterior mean of linear regression coefficients
    βhat = GeoRDD.postmean_β(mgpcv)
    # extract residuals
    residuals_rd = copy(geordd)
    residuals_rd[true]  = GeoRDD.residuals_data(geordd[true],  βhat) # treatment region
    residuals_rd[false] = GeoRDD.residuals_data(geordd[false], βhat) # control region
    # Gaussian processes fitted to residuals:
    resid_GP_dict = GeoRDD.GPRealisations(residuals_rd, mgpcv.kernel, mgpcv.logNoise)
    # obtain posterior treatment effect along border
    T,C = true, false
    sentinels = GeoRDD.sentinels(border_geo, 100)
    μpost, Σpost = GeoRDD.cliff_face(resid_GP_dict[T],
                                     resid_GP_dict[C],
                                     sentinels)
    τ_inv=GeoRDD.inverse_variance(μpost, Σpost)
    # check the credible interval is not ridiculous
    @test τ > quantile(τ_inv, 0.001)
    @test τ < quantile(τ_inv, 0.999)
    kernspace = mgpcv.kernel.kleft
    ℓ = sqrt(kernspace.ℓ2)
    maxdist = 2ℓ
    τ_proj = GeoRDD.proj_estimator(resid_GP_dict[T],
                                   resid_GP_dict[C],
                                   border_geo,
                                   2ℓ)
    @test τ > quantile(τ_proj, 0.001)
    @test τ < quantile(τ_proj, 0.999)

    # significance test
    pval_invvar_calib = GeoRDD.pval_invvar_calib(
        resid_GP_dict[T],
        resid_GP_dict[C],
        sentinels
    )
    @test pval_invvar_calib < 0.05 # effect should be detectable
end

@testset "run simulation" begin
    run_simulation()
end


text:
#!/usr/bin/env julia

#=
sum_of_odds.jl - calculate sum of odd integers
=#

__version__ = 1.1
__date__ = "2021-03-07"

up_to = 170

function is_odd(x::Int64)
  remainder = x % 2
  if remainder == 0
    return false
  else
    return true
  end
end

sum = number_of_odds = last_odd = 0

for i = 1:up_to
  if is_odd(i)
    global sum += i
    global number_of_odds += 1
    global last_odd = i
  end
end

if sum == number_of_odds ^ 2
  if number_of_odds == 0
    ;
  elseif number_of_odds < 7
    println(sum)
  else
    println("1 + 3 + 5 + ... + ", last_odd - 1,
            " + ", last_odd , " = ", sum)
  end
else
  println(stderr, "Unexpected error: ", sum)
end

max_stars_repo_name: glires/Eicosalingual

max_issues_repo_name: glires/Eicosalingual

max_forks_repo_name: glires/Eicosalingual


text:
@testset "reshaped.jl" begin
    using Distributions
    using Distributions: ReshapedDistribution
    using Test, Random, LinearAlgebra

    rng = MersenneTwister(1234)

    function test_reshaped(rng, d1, sizes)
        x1 = rand(rng, d1)
        x3 = similar(x1, size(x1)..., 3)
        rand!(rng, d1, x3)
        d1s = map(s -> reshape(d1, s...), sizes)

        # check types
        for (d, s) in zip(d1s, sizes)
            @test d isa ReshapedDistribution{length(s)}
        end

        # incorrect dimensions
        @test_throws ArgumentError reshape(d1, length(d1), 2)
        @test_throws ArgumentError reshape(d1, length(d1), 1, 2)
        @test_throws ArgumentError reshape(d1, -length(d1), -1)

        # size
        for (d, s) in zip(d1s, sizes)
            @test size(d) == s
        end

        # length
        for d in d1s
            @test length(d) == length(d1)
        end

        # rank definition for matrix distributions
        for (d, s) in zip(d1s, sizes)
            if length(s) == 2
                @test rank(d) == minimum(s)
            end
        end

        # support
        for d in d1s, s in sizes
            @test (size(d) == s) ⊻ (length(s) != length(size(d)) || !insupport(d, reshape(x1, s)))
        end

        # mean
        for (d, s) in zip(d1s, sizes)
            @test mean(d) == reshape(mean(d1), s)
        end

        # mode
        for (d, s) in zip(d1s, sizes)
            @test mode(d) == reshape(mode(d1), s)
        end

        # covariance
        for (d, s) in zip(d1s, sizes)
            @test cov(d) == cov(d1)
            if length(s) == 2
                @test cov(d, Val(false)) == reshape(cov(d1), s..., s...)
            end
        end

        # variance
        for (d, s) in zip(d1s, sizes)
            @test var(d) == reshape(var(d1), s)
        end

        # params
        for (d, s) in zip(d1s, sizes)
            @test params(d) == (d1, s)
        end

        # partype
        for d in d1s
            @test partype(d) === partype(d1)
        end

        # eltype
        for d in d1s
            @test eltype(d) === eltype(d1)
        end

        # logpdf
        for (d, s) in zip(d1s, sizes)
            @test logpdf(d, reshape(x1, s)) == logpdf(d1, x1)
            @test logpdf(d, reshape(x3, s..., 3)) == logpdf(d1, x3)
        end

        # loglikelihood
        for (d, s) in zip(d1s, sizes)
            @test loglikelihood(d, reshape(x1, s)) == loglikelihood(d1, x1)
            @test loglikelihood(d, reshape(x3, s..., 3)) == loglikelihood(d1, x3)
        end

        # rand
        for d in d1s
            x = rand(rng, d)
            @test insupport(d, x)
            @test insupport(d1, vec(x))
            @test logpdf(d, x) == logpdf(d1, vec(x))
        end

        # reshape
        for d in d1s
            @test reshape(d, size(d1)...) === d1
            @test reshape(d, size(d1)) === d1
            if d1 isa MultivariateDistribution
                @test vec(d) === d1
            end
        end
        if d1 isa MultivariateDistribution
            @test reshape(d1, size(d1)...) === d1
            @test reshape(d1, size(d1)) === d1
            @test vec(d1) === d1
        end
    end

    @testset "reshape MvNormal" begin
        σ = rand(rng, 16, 16)
        μ = rand(rng, 16)
        d1 = MvNormal(μ, σ * σ')
        sizes = [(4, 4), (8, 2), (2, 8), (1, 16), (16, 1), (4, 2, 2), (2, 4, 2), (2, 2, 2, 2)]
        test_reshaped(rng, d1, sizes)
    end

    @testset "reshape Dirichlet" begin
        α = rand(rng, 36) .+ 1 # mode is only defined if all alpha > 1
        d1 = Dirichlet(α)
        sizes = [
        (6, 6), (4, 9), (9, 4), (3, 12), (12, 3), (1, 36), (36, 1), (6, 3, 2),
        (3, 2, 6), (2, 3, 3, 2),
        ]
        test_reshaped(rng, d1, sizes)
    end

    @testset "special cases" begin
        # MatrixNormal
        rand_posdef_mat(X) = X * X' + I
        U = rand_posdef_mat(rand(5, 5))
        V = rand_posdef_mat(rand(4, 4))
        M = randn(5, 4)
        d = MatrixNormal(M, U, V)

        for v in (vec(d), reshape(d, length(d)), reshape(d, (length(d),)))
            @test v isa MvNormal
            @test mean(v) == vec(M)
            @test cov(v) == kron(V, U)
        end
    end
end


text:
"""
check_initial_conditions(initial_conditions, params)

Check whether initial conditions have been specified. Default them to standard values if not
"""
function check_initial_conditions(initial_conditions, params, grid)
if all(isnan.(initial_conditions.initial_thickness))
    default_thickness = params.default_thickness
    #@info "Did not find a specified initial thickness, reverting to default value specified in params ($default_thickness m everywhere)...\n...If you have set niter0 > 0 without invoking the update flag, you can ignore this message"
    initial_conditions = @set initial_conditions.initial_thickness =  default_thickness*ones(grid.nx, grid.ny)
end

if all(isnan.(initial_conditions.initial_viscosity))
    default_viscosity = params.default_viscosity
    #@info "Did not find a specified initial viscosity, reverting to default value specified in params ($default_viscosity Pa s everywhere)...\n...If you have set niter0 > 0 without invoking the update flag, you can ignore this message"
    initial_conditions = @set initial_conditions.initial_viscosity =  default_viscosity*ones(grid.nx, grid.ny, grid.nσ)
end

if all(isnan.(initial_conditions.initial_temperature))
    default_temperature = params.default_temperature
    #@info "Did not find a specified initial temperature, reverting to default value specified in params ($default_temperature K everywhere)...\n...If you have set niter0 > 0 without invoking the update flag, you can ignore this message"
    initial_conditions = @set initial_conditions.initial_temperature =  default_temperature*ones(grid.nx, grid.ny, grid.nσ)
end

if all(isnan.(initial_conditions.initial_damage))
    default_damage = params.default_damage
    #@info "Did not find a specified initial damage field, reverting to default value specified in params ($default_damage everywhere)...\n...If you have set niter0 > 0 without invoking the update flag, you can ignore this message"
    initial_conditions = @set initial_conditions.initial_damage =  default_damage*ones(grid.nx, grid.ny, grid.nσ)
end

#check sizes are compatible
(size(initial_conditions.initial_thickness) == (grid.nx, grid.ny)) || throw(DimensionMismatch("Initial thickness field is not compatible with grid size. Input thickess field is has size $(size(initial_conditions.initial_thickness)), which must match horizontal grid size ($(grid.nx) x $(grid.ny))"))
(size(initial_conditions.initial_temperature) == (grid.nx, grid.ny, grid.nσ)) || throw(DimensionMismatch("Initial temperature field is not compatible with grid size. Input temperature field is has size $(size(initial_conditions.initial_temperature)), which must match 3D grid size ($(grid.nx), $(grid.ny), $(grid.nσ))"))
(size(initial_conditions.initial_viscosity) == (grid.nx, grid.ny, grid.nσ)) || throw(DimensionMismatch("Initial viscosity field is not compatible with grid size. Input temperature field is has size $(size(initial_conditions.initial_temperature)), which must match 3D grid size ($(grid.nx), $(grid.ny), $(grid.nσ))"))
(size(initial_conditions.initial_damage) == (grid.nx, grid.ny, grid.nσ)) || throw(DimensionMismatch("Initial damage field is not compatible with grid size.Input temperature field is has size $(size(initial_conditions.initial_temperature)), which must match 3D grid size ($(grid.nx), $(grid.ny), $(grid.nσ))"))

return initial_conditions
end

# Utility functions
function get_bed_elevation(bed_elevation::F, grid) where (F <: Function)
bed_array = bed_elevation.(grid.xxh, grid.yyh)
return bed_array
end

function get_bed_elevation(bed_elevation::Array{T,2}, grid) where (T <: Real)
bed_array = bed_elevation
return bed_array
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
using Random, OrdinaryDiffEq

# Run this file to integrate PDE Brusselator and save snapshots.
# Configure parameters HERE, including the `saveat` parameter
# (when to save the output fields). You can choose to call this file with a
# global existing variable `u0` (the 2xNxN array of the initial condition).
# if this is the case, change the `ic` parameter below to `:given`
a = 3.0
b = 0.2
ε = 0.01
L = 300.0 # size of domain
N = 600  # no. of grid points
d = 1.0
u0 = :spiralwave
saveat = 1000:0.25:1025

# Remaining file (DO NOT ALTER!)
config = @strdict(a, b, d, L, N, saveat, u0)

function FitzHugh_Nagumo_ODE(du,u,p,t)
    a, b, d, hsq6 = p

    uu = @view u[1,:,:]
    vv = @view u[2,:,:]

    # TODO: What is this crazyiness? Creating `uuu` costs as much as computing
    # the trajectories of the solar system. But why though?
    # Aren't the Newman boundary conditions literally setting du[1, 1, 1] = du[1, N, N] = 0?

    # Create padded u matrix to incorporate Newman boundary conditions - 9-point stencil
    uuu=[ [uu[2,2] uu[2,:]' uu[2,end-1] ] ; [ uu[:,2] uu uu[:,end-1] ] ; [ uu[end-1,2] uu[end-1,:]' uu[end-1,end-1] ] ]
    diff_term = d .* (
                4 .* (uuu[2:end-1,1:end-2] .+ uuu[2:end-1,3:end] .+
                        uuu[3:end,2:end-1] .+ uuu[1:end-2,2:end-1] ) .+
                        uuu[3:end,3:end] .+ uuu[3:end,1:end-2] .+ uuu[1:end-2,3:end] .+
                        uuu[1:end-2,1:end-2] .- 20 .*uu  )  ./ hsq6

    # finite difference formula
    @. du[1,:,:] = a*uu*(1 - uu)*(uu - b) - vv + diff_term
    @. du[2,:,:] = ε*(uu - vv)
    return
end

# This is for progress tracking of DiffEq
using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

function produce_brusselator(config)
    @unpack a, b, ε, d, saveat, N, L, u0 = config
    h = L/N  # spatial step size
    hsq6 = 6*h*h

    # initial condition
    if u0 == :spiralwave
        uu = fill(0.5, N, N)
        vv = fill(0.5, N, N)
        uu[1:end, 1:(2N÷3)] .= 0
        vv[1:(2N÷3 - N÷12), 1:end] .= 0
        u0 = zeros(2,N,N)
        u0[1,:,:] = uu
        u0[2,:,:] = vv
    end
    @assert u0 isa AbstractArray
    # Solve ODE
    p = a, b, d, epsilon, hsq6
    tspan = (0.0, maximum(saveat))
    @info "Initializing"
    prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p)
    @info "Starting solve..."
    sol = solve(prob, Tsit5();
        reltol=1e-6, abstol=1e-6, maxiters = Inf, saveat,
        progress = true, progress_steps = 100,
    )
    t = sol.t
    u = sol.u
    output = @strdict u t a b d L N u0
    return output
end

prefix = length(saveat) > 10 ? "fitzhugh_dense" : "fitzhugh"

produce_or_load(datadir("FitzHugh"), config, produce_brusselator;
    prefix, suffix = "jld2", tag = false, loadfile = false,
    force = true,
)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
# Unit Test for constraint type

using COSMO, Test, SparseArrays, Random

rng = Random.MersenneTwister(1872381)

A_int = 4
b_int = 2

A_Float = 4.
b_Float = 2.

A_UInt = UInt64(4)
b_UInt = UInt64(2)

A_sparse = sprand(rng, 10, 2, 0.4)
b_sparse = sparse(rand(rng, 10))

A_rvec = [1 2 3 4]
b_rvec = 1

A_colvec = [1; 2; 3; 4]
b_colvec = [4; 3; 2; 1]

A_mat = rand(rng, 10, 10)
b_mat = sparse(rand(rng, 10, 1))

@testset "Constraints" begin

  @testset "Constructors" begin
    @test typeof(COSMO.Constraint(A_int, b_int, COSMO.ZeroSet)) <: COSMO.Constraint
    @test typeof(COSMO.Constraint(A_Float, b_Float, COSMO.ZeroSet)) <: COSMO.Constraint
    @test typeof(COSMO.Constraint(A_UInt, b_UInt, COSMO.ZeroSet)) <: COSMO.Constraint
    @test typeof(COSMO.Constraint(A_sparse, b_sparse, COSMO.ZeroSet)) <: COSMO.Constraint
    @test typeof(COSMO.Constraint(A_rvec, b_rvec, COSMO.ZeroSet)) <: COSMO.Constraint
    @test typeof(COSMO.Constraint(A_colvec, b_colvec, COSMO.ZeroSet)) <: COSMO.Constraint
    @test typeof(COSMO.Constraint(A_mat, b_mat, COSMO.ZeroSet)) <: COSMO.Constraint
  end

  @testset "Indices" begin
    A = rand(rng, 3, 3)
    b = rand(rng, 3)
    ind = 3:5
    dim = 10
    cs = COSMO.Constraint(A, b, COSMO.ZeroSet, dim, ind)

    @test cs.A[:, ind] == A
   end

  @testset "Merge Constraints" begin
    A1 = rand(rng, 10, 10)
    b1 = rand(rng, 10)
    A2 = rand(rng, 10, 10)
    b2 = rand(rng, 10)
    Am = [A1; A2]
    bm = [b1; b2]
    c1 = COSMO.Constraint(A1, b1, COSMO.ZeroSet)
    c2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)
    cArr = [c1; c2]
    cm = COSMO.Constraint(Am, bm, COSMO.ZeroSet)
    COSMO.merge_constraints!(cArr)

    @test cArr[1].A == cm.A
    @test cArr[1].b == cm.b
    @test cArr[1].convex_set.dim == cm.convex_set.dim
    @test typeof(cArr[1].convex_set) == typeof(cm.convex_set)

    c1 = COSMO.Constraint(A1, b1, COSMO.Nonnegatives)
    c2 = COSMO.Constraint(A2, b2, COSMO.Nonnegatives)
    cArr = [c1; c2]
    cm = COSMO.Constraint(Am, bm, COSMO.Nonnegatives)
    COSMO.merge_constraints!(cArr)

    @test cArr[1].A == cm.A
    @test cArr[1].b == cm.b
    @test cArr[1].convex_set.dim == cm.convex_set.dim
    @test typeof(cArr[1].convex_set) == typeof(cm.convex_set)
  end

end
nothing


text:
# Buffer
# ======

# Data Layout
# -----------
#
# Buffered data are stored in `data` and three position fields are used to keep
# track of marked data, buffered data and margin.
#
#                     marked      buffer      margin
#                  |<-------->||<-------->||<-------->|
#       data   ....xxxxxxxxxxxxXXXXXXXXXXXX............
#              ^   ^           ^           ^          ^
#   position   1   markpos     bufferpos   marginpos  lastindex(data)
#
# `markpos` is positive iff there are marked data; otherwise it is set to zero.
# `markpos` ≤ `bufferpos` ≤ `marginpos` must hold whenever possible.

mutable struct Buffer
    # data and positions (see above)
    data::Vector{UInt8}
    markpos::Int
    bufferpos::Int
    marginpos::Int

    # the total number of transcoded bytes
    transcoded::Int64

    function Buffer(size::Integer)
        return new(Vector{UInt8}(undef, size), 0, 1, 1, 0)
    end

    function Buffer(data::Vector{UInt8})
        return new(data, 0, 1, length(data)+1, 0)
    end
end

function Buffer(data::Base.CodeUnits{UInt8})
    return Buffer(Vector{UInt8}(data))
end

function Base.length(buf::Buffer)
    return length(buf.data)
end

function bufferptr(buf::Buffer)
    return pointer(buf.data, buf.bufferpos)
end

function buffersize(buf::Buffer)
    return buf.marginpos - buf.bufferpos
end

function buffermem(buf::Buffer)
    return Memory(bufferptr(buf), buffersize(buf))
end

function marginptr(buf::Buffer)
    return pointer(buf.data, buf.marginpos)
end

function marginsize(buf::Buffer)
    return lastindex(buf.data) - buf.marginpos + 1
end

function marginmem(buf::Buffer)
    return Memory(marginptr(buf), marginsize(buf))
end

function ismarked(buf::Buffer)
    return buf.markpos != 0
end

function mark!(buf::Buffer)
    return buf.markpos = buf.bufferpos
end

function unmark!(buf::Buffer)
    if buf.markpos == 0
        return false
    else
        buf.markpos = 0
        return true
    end
end

function reset!(buf::Buffer)
    @assert buf.markpos > 0
    buf.bufferpos = buf.markpos
    buf.markpos = 0
    return buf.bufferpos
end

# Notify that `n` bytes are consumed from `buf`.
function consumed!(buf::Buffer, n::Integer; transcode::Bool = false)
    buf.bufferpos += n
    if transcode
        buf.transcoded += n
    end
    return buf
end

# Notify that `n` bytes are supplied to `buf`.
function supplied!(buf::Buffer, n::Integer; transcode::Bool = false)
    buf.marginpos += n
    if transcode
        buf.transcoded += n
    end
    return buf
end

# Discard buffered data and initialize positions.
function initbuffer!(buf::Buffer)
    buf.markpos = buf.transcoded = 0
    buf.bufferpos = buf.marginpos = 1
    return buf
end

# Remove all buffered data.
function emptybuffer!(buf::Buffer)
    buf.marginpos = buf.bufferpos
    return buf
end

# Make margin with ≥`minsize` and return the size of it.
# If eager is true, it tries to move data even when the buffer has enough margin.
function makemargin!(buf::Buffer, minsize::Integer; eager::Bool = false)
    @assert minsize ≥ 0
    if buffersize(buf) == 0 && buf.markpos == 0
        buf.bufferpos = buf.marginpos = 1
    end
    if marginsize(buf) < minsize || eager
        # shift data to left
        if buf.markpos == 0
            datapos = buf.bufferpos
            datasize = buffersize(buf)
        else
            datapos = buf.markpos
            datasize = buf.marginpos - buf.markpos
        end
        copyto!(buf.data, 1, buf.data, datapos, datasize)
        shift = datapos - 1
        if buf.markpos > 0
            buf.markpos -= shift
        end
        buf.bufferpos -= shift
        buf.marginpos -= shift
    end
    if marginsize(buf) < minsize
        # expand data buffer
        resize!(buf.data, buf.marginpos + minsize - 1)
    end
    @assert marginsize(buf) ≥ minsize
    return marginsize(buf)
end

# Read a byte.
function readbyte!(buf::Buffer)
    b = buf.data[buf.bufferpos]
    consumed!(buf, 1)
    return b
end

# Write a byte.
function writebyte!(buf::Buffer, b::UInt8)
    buf.data[buf.marginpos] = b
    supplied!(buf, 1)
    return 1
end

# Skip `n` bytes in the buffer.
function skipbuffer!(buf::Buffer, n::Integer)
    @assert n ≤ buffersize(buf)
    consumed!(buf, n)
    return buf
end

# Take the ownership of the marked data.
function takemarked!(buf::Buffer)
    @assert buf.markpos > 0
    sz = buf.marginpos - buf.markpos
    copyto!(buf.data, 1, buf.data, buf.markpos, sz)
    initbuffer!(buf)
    return resize!(buf.data, sz)
end

# Copy data from `data` to `buf`.
function copydata!(buf::Buffer, data::Ptr{UInt8}, nbytes::Integer)
    makemargin!(buf, nbytes)
    unsafe_copyto!(marginptr(buf), data, nbytes)
    supplied!(buf, nbytes)
    return buf
end

# Copy data from `buf` to `data`.
function copydata!(data::Ptr{UInt8}, buf::Buffer, nbytes::Integer)
    # NOTE: It's caller's responsibility to ensure that the buffer has at least
    # nbytes.
    @assert buffersize(buf) ≥ nbytes
    unsafe_copyto!(data, bufferptr(buf), nbytes)
    consumed!(buf, nbytes)
    return data
end

# Insert data to the current buffer.
function insertdata!(buf::Buffer, data::Ptr{UInt8}, nbytes::Integer)
    makemargin!(buf, nbytes)
    copyto!(buf.data, buf.bufferpos + nbytes, buf.data, buf.bufferpos, buffersize(buf))
    unsafe_copyto!(bufferptr(buf), data, nbytes)
    supplied!(buf, nbytes)
    return buf
end

# Find the first occurrence of a specific byte.
function findbyte(buf::Buffer, byte::UInt8)
    p = ccall(
        :memchr,
        Ptr{UInt8},
        (Ptr{UInt8}, Cint, Csize_t),
        pointer(buf.data, buf.bufferpos), byte, buffersize(buf))
    if p == C_NULL
        return marginptr(buf)
    else
        return p
    end
end


text:
# timing

# system date in seconds
time() = ccall(:clock_now, Float64, ())

# high-resolution relative time, in nanoseconds
time_ns() = ccall(:jl_hrtime, Uint64, ())

# total number of bytes allocated so far
gc_bytes() = ccall(:jl_gc_total_bytes, Int64, ())

function tic()
    t0 = time_ns()
    task_local_storage(:TIMERS, (t0, get(task_local_storage(), :TIMERS, ())))
    return t0
end

function toq()
    t1 = time_ns()
    timers = get(task_local_storage(), :TIMERS, ())
    if is(timers,())
        error("toc() without tic()")
    end
    t0 = timers[1]::Uint64
    task_local_storage(:TIMERS, timers[2])
    (t1-t0)/1e9
end

function toc()
    t = toq()
    println("elapsed time: ", t, " seconds")
    return t
end

# print elapsed time, return expression value
macro time(ex)
    quote
        local b0 = gc_bytes()
        local t0 = time_ns()
        local val = $(esc(ex))
        local t1 = time_ns()
        local b1 = gc_bytes()
        println("elapsed time: ", (t1-t0)/1e9, " seconds (", b1-b0, " bytes allocated)")
        val
    end
end

# print nothing, return elapsed time
macro elapsed(ex)
    quote
        local t0 = time_ns()
        local val = $(esc(ex))
        (time_ns()-t0)/1e9
    end
end

# measure bytes allocated without any contamination from compilation
macro allocated(ex)
    quote
        let
            local f
            function f()
                b0 = gc_bytes()
                $(esc(ex))
                b1 = gc_bytes()
                b1-b0
            end
            f()
        end
    end
end

# print nothing, return value, elapsed time & bytes allocated
macro timed(ex)
    quote
        local b0 = gc_bytes()
        local t0 = time_ns()
        local val = $(esc(ex))
        local t1 = time_ns()
        local b1 = gc_bytes()
        val, (t1-t0)/1e9, b1-b0
    end
end

# BLAS utility routines
function blas_vendor()
    try
        cglobal((:openblas_set_num_threads, Base.libblas_name), Void)
        return :openblas
    end
    try
        cglobal((:MKL_Set_Num_Threads, Base.libblas_name), Void)
        return :mkl
    end
    return :unknown
end

openblas_get_config() = chop(bytestring( ccall((:openblas_get_config, Base.libblas_name), Ptr{Uint8}, () )))

function blas_set_num_threads(n::Integer)
    blas = blas_vendor()
    if blas == :openblas
        return ccall((:openblas_set_num_threads, Base.libblas_name), Void, (Int32,), n)
    elseif blas == :mkl
        # MKL may let us set the number of threads in several ways
        return ccall((:MKL_Set_Num_Threads, Base.libblas_name), Void, (Cint,), n)
    end

    # OSX BLAS looks at an environment variable
    @osx_only ENV["VECLIB_MAXIMUM_THREADS"] = n

    return nothing
end

function check_blas()
    blas = blas_vendor()
    if blas == :openblas
        openblas_config = openblas_get_config()
        openblas64 = ismatch(r".*USE64BITINT.*", openblas_config)
        if Base.USE_BLAS64 != openblas64
            if !openblas64
                println("ERROR: OpenBLAS was not built with 64bit integer support.")
                println("You're seeing this error because Julia was built with USE_BLAS64=1")
                println("Please rebuild Julia with USE_BLAS64=0")
            else
                println("ERROR: Julia was not built with support for OpenBLAS with 64bit integer support")
                println("You're seeing this error because Julia was built with USE_BLAS64=0")
                println("Please rebuild Julia with USE_BLAS64=1")
            end
            println("Quitting.")
            quit()
        end
    elseif blas == :mkl
        if Base.USE_BLAS64
            ENV["MKL_INTERFACE_LAYER"] = "ILP64"
        end
    end

    #
    # Check if BlasInt is the expected bitsize, by triggering an error
    #
    (_, info) = LinAlg.LAPACK.potrf!('U', [1.0 0.0; 0.0 -1.0])
    if info != 2 # mangled info code
        if info == 2^33
            error("""BLAS and LAPACK are compiled with 32-bit integer support, but Julia expects 64-bit integers. Please build Julia with USE_BLAS64=0.""")
        elseif info == 0
            error("""BLAS and LAPACK are compiled with 64-bit integer support but Julia expects 32-bit integers. Please build Julia with USE_BLAS64=1.""")
        else
            error("""The LAPACK library produced an undefined error code. Please verify the installation of BLAS and LAPACK.""")
        end
    end

end

## printing with color ##

function with_output_color(f::Function, color::Symbol, io::IO, args...)
    buf = IOBuffer()
    have_color && print(buf, get(text_colors, color, color_normal))
    try f(buf, args...)
    finally
        have_color && print(buf, color_normal)
        print(io, takebuf_string(buf))
    end
end

print_with_color(color::Symbol, io::IO, msg::String...) =
    with_output_color(print, color, io, msg...)
print_with_color(color::Symbol, msg::String...) =
    print_with_color(color, STDOUT, msg...)
println_with_color(color::Symbol, io::IO, msg::String...) =
    with_output_color(println, color, io, msg...)
println_with_color(color::Symbol, msg::String...) =
    println_with_color(color, STDOUT, msg...)

## warnings and messages ##

function info(msg::String...; prefix="INFO: ")
    println_with_color(:blue, STDERR, prefix, chomp(string(msg...)))
end

# print a warning only once

const have_warned = Set()
warn_once(msg::String...) = warn(msg..., once=true)

function warn(msg::String...; prefix="WARNING: ", once=false, key=nothing, bt=nothing)
    str = chomp(bytestring(msg...))
    if once
        if key === nothing
            key = str
        end
        (key in have_warned) && return
        push!(have_warned, key)
    end
    with_output_color(:red, STDERR) do io
        print(io, prefix, str)
        if bt !== nothing
            show_backtrace(io, bt)
        end
        println(io)
    end
end

warn(err::Exception; prefix="ERROR: ", kw...) =
    warn(sprint(io->showerror(io,err)), prefix=prefix; kw...)

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
using Documenter, SossModels

makedocs(;
    modules=[SossModels],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/cscherrer/SossModels.jl/blob/{commit}{path}#L{line}",
    sitename="SossModels.jl",
    authors="Chad Scherrer",
    assets=String[],
)

deploydocs(;
    repo="github.com/cscherrer/SossModels.jl",
)


text:
using JML
using Base.Test

# write your own tests here

anyerrors = false

testcases = [
            "FunctionsTest.jl"
           ]

Base.info("Running tests...")

for case in testcases
    try
        include(case)
        println("\t\033[1m\033[32mPASSED\033[0m: $(case)")
    catch e
        anyerrors = true
        println("\t\033[1m\033[31mFAILED\033[0m: $(case)")
    end
end

if anyerrors
    throw("Test failed")
end


text:
export fit!, project, conmatrix

#=
Maximum Liklihood

Linear Discriminant Analysis
Quadratic Discriminant Analysis
Diagonal Discriminant Analysis
Diagonal Quadratic Discriminant Analysis
Shrinkage?
=#

#=
LDA
=#

function fit!{C<:LDA, V<:validation}(m::decoder{C,V},response::Array{Float64,2},stimulus::Array{Float64,1})

    m.classes=unique(stimulus)
    k=length(m.classes)

    #if any stimulus appears <10 times, ignore it
    keep=Array(Int64,0)
    for i=1:k
        Group = sum(stimulus.==m.classes[i])
        if Group<10
        else
            push!(keep,i)
        end
    end

    m.classes=m.classes[keep]
    k=length(m.classes)

    nGroup=zeros(Int64,k)
    GroupMean=zeros(Float64,k,size(response,2))
    Sw=zeros(Float64,size(response,2),size(response,2))
    m.c.W=zeros(Float64,k,size(response,2))
    covmat=zeros(Float64,size(response,2),size(response,2),k)

    for i=1:k
        Group = stimulus.==m.classes[i]
        nGroup[i]=sum(Group)

        GroupMean[i,:]=mean(response[Group,:],1)

        covmat[:,:,i] = cov(response[Group,:])

    end

    for i=1:k
        Sw+=(nGroup[i]/sum(nGroup)).*covmat[:,:,i]
    end

    St=cov(response)
    Sb = St - Sw

    (m.c.myv, W)=eig(Sb,Sw)

    mostvar=find((m.c.myv./sum(m.c.myv)).>.01)

    m.c.W=W[:,mostvar]

    m.c.centroids=GroupMean*m.c.W

    nothing

end

function project{C<:LDA, V<:validation}(m::decoder{C,V},response::Array{Float64,2},stimulus::Array{Float64,1})

    predict=zeros(Float64,size(stimulus,1))

    xnew=response*m.c.W

    mydist=zeros(Float64,length(m.classes))

    for i=1:size(stimulus,1)
        for j=1:length(m.classes)
            mydist[j]=norm(xnew[i,:]-m.c.centroids[j,:])
        end
        predict[i]=m.classes[indmin(mydist)]
    end

    predict

end

function project(myLDA::PyObject,response::Array{Float64,2},stimulus::Array{Float64,1})

    predict=zeros(Float64,size(stimulus,1))

    classes=myLDA[:classes_]
    mydist=zeros(Float64,length(classes))
    mycentroids=myLDA[:transform](myLDA[:means_])

    for i=1:size(stimulus,1)
        xnew=myLDA[:transform](response[i,:])
        for j=1:length(classes)
            mydist[j]=norm(xnew-mycentroids[j,:])
        end
        predict[i]=classes[indmin(mydist)]
    end

    predict

end

function validate{C<:LDA, V<:Training}(m::decoder{C,V},response::Array{Float64,2},stimulus::Array{Float64,1})

end

function validate{C<:LDA, V<:LeaveOne}(m::decoder{C,V},response::Array{Float64,2},stimulus::Array{Float64,1})

end

#=
Nearest Neighbor Classification

Quian Quiroga et al 2006
=#

#=
Naive Bayesian
=#

#=
General Methods
=#

function conmatrix{C<:classifier,V<:validation}(m::decoder{C,V},predict::Array{Float64,1},stimulus::Array{Float64,1})

    inds=collect(1:length(m.classes))

    mydict=Dict{Float64,Int64}([m.classes[i]=>inds[i] for i=1:length(m.classes)])

    #y is true, x is predicted
    #sum across columns should equal 1
    conmat=zeros(Float64,length(m.classes),length(m.classes))

    for i=1:length(stimulus)
        if length(find(stimulus[i].==m.classes))==0 #might be a stimulus in a val set that isn't in training
        else
            yind=mydict[stimulus[i]]
            xind=mydict[predict[i]]
            conmat[yind,xind]+=1
        end
    end

    totals_real=sum(conmat,2)
    totals_predict=sum(conmat,1)

    for i=1:length(m.classes)
        conmat[i,:]=conmat[i,:]./totals_real[i]
    end

    p_real=squeeze(totals_real./sum(totals_real),2)
    p_predict=squeeze((totals_predict./sum(totals_predict))',2)

    (conmat, p_real, p_predict)

end


text:
@testset "HW4" begin
    using DMUStudent.HW4: gw
    using CommonRLInterface
    using Compose

    A = actions(gw)
    done = false
    rsum = 0.0
    reset!(gw)
    s = observe(gw)
    t = 1
    γ = 0.95
    max_t = 200
    while t <= max_t && !terminated(gw)
        a = rand(A)

        r = act!(gw, a)
        sp = observe(gw)

        rsum += γ^t*r

        s = sp
        t += 1
    end
    @test -12.0 <= rsum <= 10.0

    # rendering gw
    rndr = HW4.render(gw, color=s->10.0*rand(), policy=s->rand(A))
    @test rndr isa Context
    sprint((io, x)->show(io, MIME("text/html"), x), rndr)

    rndr = HW4.render(gw)
    @test rndr isa Context
    sprint((io, x)->show(io, MIME("text/html"), x), rndr)
end