text:
# trait for iterablility
abstract type IterableTrait end
struct Iterable{T} <: IterableTrait end
struct NotIterable{T} <: IterableTrait end
IterableTrait(::Type) = NotIterable{Any}()
IterableTrait(::Type{A}) where {A<:AbstractArray{T}} where {T} = Iterable{T}()
IterableTrait(::Type{A}) where {A<:Dict{U,T}} where {T} where {U} = Iterable{T}()

# definition of a tree
abstract type AbstractTree end

mutable struct Leaf <: AbstractTree
    ID::Union{Nothing,AbstractTree}
    name::String
    parent::Union{Nothing,AbstractTree}
    ext::Dict{Symbol,Any}
    function Leaf()
        return new(nothing,"",nothing,Dict{Symbol,Any}())
    end
    function Leaf(parent::AbstractTree)
        return new(nothing,"",parent,Dict{Symbol,Any}())
    end
    function Leaf(name::String)
        return new(nothing,name,nothing,Dict{Symbol,Any}())
    end
    function Leaf(name::String,parent::AbstractTree)
        return new(nothing,name,parent,Dict{Symbol,Any}())
    end
end

mutable struct Tree <: AbstractTree
    children::Union{Array{AbstractTree,1},Array{Leaf,1},Array{Tree,1}}
    ID::Union{Nothing,AbstractTree}
    name::String
    parent::Union{Nothing,AbstractTree}
    ext::Dict{Symbol,Any}
    Tree(children) = Tree(IterableTrait(typeof(children)), children)
    function Tree(::Iterable{<:AbstractTree}, children)
        new(children, nothing, "",nothing,Dict{Symbol,Any}())
    end
    function Tree()
        new(Array{AbstractTree,1}(),nothing,"",nothing,Dict{Symbol,Any}())
    end
    function Tree(parent::AbstractTree)
        new(Array{AbstractTree,1}(),nothing,"",parent,Dict{Symbol,Any}())
    end
    function Tree(name::String)
        new(Array{AbstractTree,1}(),nothing,name,nothing,Dict{Symbol,Any}())
    end
    function Tree(name::String,parent::AbstractTree)
        new(Array{AbstractTree,1}(),nothing,name,parent,Dict{Symbol,Any}())
    end
end

function Base.copy(leaf::Leaf)
    Leaf()
end

function Base.copy(tree::Tree)
    Tree(map(copy, tree.children))
end

function Base.map(f, leaf::Leaf)
    f(leaf)
end

function Base.map(f, tree::Tree)
    Tree(map(x -> map(f, x), tree.children))
end

function Base.map(f, dict::Dict)
    Dict(key => f(dict[key]) for key in keys(dict))
end

function Base.show(io::IO, tree::AbstractTree)
    if typeof(tree) == Tree
        if tree.name != ""
            print(
                io,
                "Subtree rooted at node " *
                tree.name *
                " containing " *
                string(count(tree)) *
                " nodes",
            )
        else
            print(io, "Subtree containing " * string(count(tree)) * " nodes")
        end
    else
        if tree.name != ""
            print(io, "Leaf node " * tree.name)
        else
            print(io, "Leaf node")
        end
    end
end

function getID(node::AbstractTree)
    if node.ID==nothing
        return node
    else
        return node.ID
    end
end

function count(tree::AbstractTree)
    i = 1
    if typeof(tree) == Tree
        for child in tree.children
            i += count(child)
        end
    end
    return i
end

function print_tree(some_tree::AbstractTree)
    function helper(tree::Tree, depth)
        println("  "^depth * "--" * tree.name)
        for child in tree.children
            helper(child, depth + 1)
        end
    end
    function helper(leaf::Leaf, depth)
        println("  "^depth * "--" * leaf.name)
    end
    helper(some_tree, 0)
    nothing
end

"""
	print_tree(some_tree::AbstractTree, data::Dict{AbstractTree,Any})

Given `some_tree`, this function prints a simple representation of the tree to the REPL.

### Required Arguments
`some_tree` is the tree we wish to visualise

### Optional Arguments
`data` is a dictionary indexed by each node in `some_tree`
"""
function print_tree(some_tree::AbstractTree, data::Dict{AbstractTree,T} where T <: Any)
    function helper(tree::Tree, depth)
        if tree in keys(data)
            println("  "^depth * "--" * tree.name * " (" * string(data[tree]) * ")")
        elseif getID(tree) in keys(data)
            println("  "^depth * "--" * tree.name * " (" * string(data[getID(tree)]) * ")")
        else
            println("  "^depth * "--" * tree.name)
        end

        for child in tree.children
            helper(child, depth + 1)
        end
    end
    function helper(leaf::Leaf, depth)
        if leaf in keys(data)
            println("  "^depth * "--" * leaf.name * " (" * string(data[leaf]) * ")")
        elseif getID(leaf) in keys(data)
            println("  "^depth * "--" * leaf.name * " (" * string(data[getID(leaf)]) * ")")
        else
            println("  "^depth * "--" * leaf.name)
        end
    end
    helper(some_tree, 0)

    nothing
end

"""
	visualize_tree(some_tree::AbstractTree,
        data::Union{Dict{Symbol,Any},Dict{Symbol,Dict{AbstractTree,Float64}}};
        scale_edges=nothing,
        scale_all=1.0)

Given `some_tree`, this function generates a html/js visualization of the tree.

### Required Arguments
`some_tree` is the tree we wish to visualise.

`data` is a dictionary of the data we wish to display, each element is another
dictionary indexed by the nodes of the tree.

### Optional Arguments
`scale_edges` this scales the lengths of the arcs.

`scale_all` this scales the whole network.
"""
function visualize_tree(some_tree::AbstractTree, data::Union{Dict{Symbol,Any},Dict{Symbol,Dict{AbstractTree,Float64}}};scale_edges=nothing,scale_all=1.0)
    maxdata=Dict{Symbol,Any}()
    mindata=Dict{Symbol,Any}()

    function node_json(node)
        temp="{"
        temp*="id:"
        temp*=string(get_id[node])
        temp*=",label:\""
        temp*=node.name*"\""
        temp*=",posX:"
        temp*=string(position[node][1])
        temp*=",posY:"
        temp*=string(position[node][2])
        temp*=",data:{"
        first=true
        for sym in keys(data)
            if typeof(collect(keys(data[sym]))[1]) <: AbstractTree
                if node in keys(data[sym])
                    if first
                        first=false
                    else
                        temp*=","
                    end
                    temp*=string(sym)
                    temp*=":"
                    if typeof(data[sym][node])==String
                        temp*="\""
                        temp*=data[sym][node]
                        temp*="\""
                    else
                        if sym ∉ keys(maxdata)
                            maxdata[sym]=data[sym][node]
                            mindata[sym]=data[sym][node]
                        else
                            if data[sym][node] < mindata[sym]
                                mindata[sym]=data[sym][node]
                            elseif data[sym][node] > maxdata[sym]
                                maxdata[sym]=data[sym][node]
                            end
                        end
                        temp*=string(data[sym][node])
                    end
                end
            elseif typeof(collect(keys(data[sym]))[1]) == String
                first2=true
                for key in keys(data[sym])
                    if node in keys(data[sym][key])
                        if first
                            first=false
                        else
                            temp*=","
                        end
                        if first2
                            first2=false
                            # temp*=","
                            temp*=string(sym)
                            temp*=":{"
                        end
                        temp*=""
                        temp*=replace(key, "," => "_")
                        temp*=":"
                        if typeof(data[sym][key][node])==String
                            temp*="\""
                            temp*=data[sym][key][node]
                            temp*="\""
                        else
                            if sym ∉ keys(maxdata)
                                maxdata[sym]=data[sym][key][node]
                                mindata[sym]=data[sym][key][node]
                            else
                                if data[sym][key][node] < mindata[sym]
                                    mindata[sym]=data[sym][key][node]
                                elseif data[sym][key][node] > maxdata[sym]
                                    maxdata[sym]=data[sym][key][node]
                                end
                            end
                            temp*=string(data[sym][key][node])
                        end
                    end
                end
                if !first2
                    temp*="}"
                end
            end
        end
        temp*="},level:"
        temp*=string(depth(node))
        temp*=",},"
    end

    function arc_json(node,parent)
        temp="{"
        temp*="from:"
        temp*=string(get_id[parent])
        temp*=",to:"
        temp*=string(get_id[node])
        temp*=",},"
    end

    index=1
    nodes="var nodes = ["
    arcs="var edges = ["
    get_id=Dict{AbstractTree,Int64}()

    angles=Dict{AbstractTree,Float64}()
    position=Dict{AbstractTree,Tuple{Float64,Float64}}()

    function setpositions(node::AbstractTree,alpha::Float64,l::Float64,scale::Float64;first=false)
        if typeof(node)==Leaf
            return
        end
        a=2*pi/(length(node.children)-1+2*alpha)

        current=angles[node]+pi+a*alpha
        for child in node.children
            angles[child]=current
            position[child]=(position[node][1]+l*cos(current),position[node][2]+l*sin(current))
            setpositions(child,alpha,l*scale,scale)
            current+=a
        end
    end

    scale_factors = [ [1.0],
                      [1.0,0.9,0.85,0.78,0.73,0.7,0.67,0.65],
                      [1.0,0.65,0.52,0.48],
                      [1.0,0.45,0.42,0.42,0.41,0.4],
                      [1.0,0.44,0.37,0.36],
                      [1.0,0.42,0.34,0.33],
                      [1.0,0.35,0.3],
                      [1.0,0.3,0.26],
                      [1.0,0.27,0.23],
                      [1.0,0.24,0.22] ]

    if scale_edges==nothing
        if typeof(some_tree)==Leaf
            scale_edges=1.0
        else
            dg=length(some_tree.children)
            if dg<=10
                dp=min(depth(collect(some_tree)[end]),length(scale_factors[dg]))
                scale_edges=scale_factors[dg][dp]
            else
                scale_edges=0.22*0.91^(dg-10)
            end
        end
    end

    angles[some_tree]=0.0
    position[some_tree]=(0.0,0.0)

    setpositions(some_tree,0.5,700.0*scale_all,scale_edges,first=true)

    for node in collect(some_tree)
        get_id[node]=index
        nodes*=node_json(node)
        parent=node.parent
        if parent!=nothing
            arcs*=arc_json(node,parent)
        end
        index+=1
    end

    nodes*="];"
    arcs*="];"

    scale="var scale = {"

    for sym in keys(maxdata)
        if length(scale)!=13
            scale*=","
        end
        scale*=string(sym)
        scale*=":"
        scale*="{min:"
        scale*=string(mindata[sym])
        scale*=",max:"
        if maxdata[sym]==mindata[sym]
            maxdata[sym]+=1
        end
        scale*=string(maxdata[sym])
        scale*="}"
    end
    scale*="};"
    min_size=(1.0+79.0/scale_all)*(scale_edges^(0.65/scale_all))^(depth(collect(some_tree)[end]))*scale_all
    max_size=(1.0+79.0/scale_all)*(scale_all)
    data=">"*nodes*"\n"*arcs*"\n"*scale*"\n"*"var node_scale="*string(scale_edges^(0.65/scale_all))*";"*"var min_size="*string(min_size)*";"*"var max_size="*string(max_size)*";"

    s = read(joinpath(dirname(@__DIR__),"visualise","visualize.html"), String)
    c = ">"*read(joinpath(dirname(@__DIR__),"visualise","colors.js"), String)
    filename=joinpath("vis"*string(Int(round((time()*100)%100000)))*".html")
    file=open(filename,"w")
    println(file,replace(replace(s," src=\"smallnetwork.js\">" => data)," src=\"colors.js\">" => c))
    close(file)

    if Sys.iswindows()
        run(`$(ENV["COMSPEC"]) /c start $(filename)`)
    elseif Sys.isapple()
        run(`open $(filename)`)
    elseif Sys.islinux() || Sys.isbsd()
        run(`xdg-open $(filename)`)
    else
        error("Unable to show plot. Try opening the file $(filename) manually.")
    end
    return
end

"""
	collect(tree::Tree;order=:depth)

Given `tree`, this function returns an array of corresponding nodes. By default this will be in a depth-first order.

### Required Arguments
`tree` is the tree from which we wish to collect the nodes

### Optional Arguments
`order` can be set to `:depth` or `:breadth` to specify the order that the nodes are listed in the array.

### Examples
    nodes = collect(tree) #gets an array of nodes from tree in depth-first order
    nodes = collect(tree,order=:breadth) #gets an array of nodes from tree in breadth-first order
"""
function Base.collect(tree::Tree;order=:depth)
    index=1
    collection=Array{AbstractTree,1}()
    push!(collection,tree)
    while index<=length(collection)
        if typeof(collection[index])==Tree
            if order==:depth
                for i in eachindex(collection[index].children)
                    insert!(collection,index+i,collection[index].children[i])
                end
            elseif order==:breadth
                append!(collection,collection[index].children)
            else
                error("\'order\' should be set to :depth or :breadth")
            end
        end
        index+=1
    end
    collection
end

function Base.collect(leaf::Leaf;order=:depth)
    [leaf]
end

function label_nodes(tree::AbstractTree)
    for t in collect(tree)
        t.name=history2(t)
    end
    tree
end

"""
	get_leafnodes(tree::AbstractTree)

Given `tree`, this function returns an array of corresponding `Leaf` nodes.

### Required Arguments
`tree` is the tree from which we wish to collect leaf nodes

### Example
    leafnodes = JuDGE.get_leafnodes(tree)
"""
function get_leafnodes(tree::AbstractTree)
    function helper(leaf::Leaf, collection)
        push!(collection, leaf)
    end
    function helper(someTree::Tree, collection)
        for child in someTree.children
            helper(child, collection)
        end
    end
    result = Array{Leaf,1}()
    helper(tree, result)
    result
end

"""
	ConditionallyUniformProbabilities(tree::AbstractTree)

Given a tree, this function returns a dictionary which maps nodes of the tree to probabilities, given that there are conditionally uniform probabilities over the children of any node.

### Required Arguments
`tree` is the tree for which the probability distribution will be generated

### Example
    probs = ConditionallyUniformProbabilities(tree)
"""
function ConditionallyUniformProbabilities(tree::T where {T<:AbstractTree})
    function result(subtree::T where {T<:AbstractTree})
        if subtree == tree
            return 1.0
        else
            parent = subtree.parent
            return result(parent) / length(parent.children)
        end
    end
    prob=Dict{AbstractTree,Float64}()
    for node in collect(tree)
        prob[node]=result(node)
    end
    return prob
end

"""
	UniformLeafProbabilities(tree::AbstractTree)

Given a tree, this function returns a dictionary which maps nodes of the tree to probabilities, given that there is a uniform distribution over the leaf nodes

### Required Arguments
`tree` is the tree for which the probability distribution will be generated

### Example
    probs = UniformLeafProbabilities(tree)
"""
function UniformLeafProbabilities(tree::T where {T<:AbstractTree})
    function result(subtree::T where {T<:AbstractTree})
        if typeof(subtree) == Leaf
            return p
        else
            pr=0
            for child in subtree.children
                pr+=result(child)
            end
            return pr
        end
    end

    leafnodes=get_leafnodes(tree)
    p=1/length(leafnodes)
    prob=Dict{AbstractTree,Float64}()
    for node in collect(tree)
        prob[node]=result(node)
    end
    return prob
end

"""
	convert_probabilities(tree::AbstractTree, probabilities::Dict{NodeID,Float64})

Given a dictionary of conditional probabilities for each node in `tree`, this function
returns a dictionary that maps the `NodeID` of each node of `tree` to the corresponding
unconditional probability.

### Required Arguments
`tree` is the tree that the probabilities pertain to

`probabilities` is a dictionary of condition probabilities for each node in `tree`

### Example
    probs = JuDGE.convert_probabilities(tree,probabilities)
"""
function convert_probabilities(tree::T where {T<:AbstractTree}, probabilities::Dict{AbstractTree,Float64})
    function result(subtree::T where {T<:AbstractTree})
        if subtree == tree
            return probabilities[getID(subtree)]
        else
            parent = subtree.parent
            return result(parent) * probabilities[getID(subtree)]
        end
    end
    prob=Dict{AbstractTree,Float64}()
    for node in collect(tree)
        prob[node]=result(node)
    end
    return prob
end

"""
	depth(tree::AbstractTree)

Given `tree`, this function returns the depth.
The root node has a depth of 0.

### Required Arguments
`tree` is the node we wish to find the depth for.

### Example
    depth = JuDGE.depth(tree) #returns the depth of a node in a tree
"""
function depth(tree::AbstractTree)
    length(history(tree))-1
end

"""
	history(tree::AbstractTree)

Given `tree`, this function returns the history back up to the root node of `tree`.

### Required Arguments
`tree` is the node that we wish to find the history for.

### Example
    history = JuDGE.history(tree) #get a vector of nodes that precede tree
"""
function history(tree::T where {T<:AbstractTree})
    function helper(state, subtree)
        if subtree.parent == nothing
            push!(state, subtree)
            return state
        else
            push!(state, subtree)
            helper(state, subtree.parent)
        end
        state
    end
    helper(Array{AbstractTree,1}(), tree)
end

function history2(tree::T where {T<:AbstractTree})
    function helper(state, subtree)
        if subtree.parent == nothing
            state = "1" * state
            return state
        else
            parent = subtree.parent
            for i = 1:length(parent.children)
                if subtree == parent.children[i]
                    state = string(i) * state
                    break
                end
            end
            state = helper(state, subtree.parent)
        end
        state
    end
    helper("", tree)
end

"""
	narytree(depth::Int64, degree::Int64)

Given the `depth` and `degree`, this function returns an N-ary tree. Note that a depth of 0 return a single `Leaf` node (which is also the root node of the tree).

### Required Arguments
`depth` is the maximum number of arcs from the root node any `Leaf` node

`degree` is the number of children of all nodes, other than the `Leaf` nodes

### Example
    tree = narytree(2,2)
"""
function narytree(depth::Int64, degree::Int64)
    if depth==0
        return Leaf("1")
    end

    tree=Tree("1")
    layer=[tree]
    for dp in 1:depth
        nextlayer=AbstractTree[]
        for t in layer
            for dg in 1:degree
                if dp!=depth
                    nt=Tree(t.name*string(dg),t)
                    push!(t.children,nt)
                    push!(nextlayer,nt)
                else
                    push!(t.children,Leaf(t.name*string(dg),t))
                end
            end
        end
        layer=nextlayer
    end
    #label_nodes(tree)
    tree
end

"""
	get_node(tree::AbstractTree, indices::Array{Int64,1})

Given a `tree`, and an array of `indices`, this function returns the corresponding node in the tree.

### Required Arguments
`tree` is the tree from which we are finding the node

`indicies` is an array of integer indices identifying a node within `tree`.

### Examples
    node = get_node(tree,[1]) #get the root node
    node = get_node(tree,[1,1]) #get the first child of the root node
    node = get_node(tree,[1,2]) #get the second child of the root node
"""
function get_node(tree::AbstractTree, indices::Array{Int64,1})
    node = tree
    for i = 2:length(indices)
        node = node.children[indices[i]]
    end
    return node
end

"""
	get_groups(tree::AbstractTree; combine=0)

Given a `tree`, this function will split it up into an array of subtrees. These can
be provided as `blocks` for the `JuDGE.solve()` function to perform partial pricing.

### Required Arguments
`tree` is the tree from which we are finding the node

### Optional Arguments
`combine` this parameter determines the size of the subtrees (higher creates larger subtrees).
If set to 0, the subtrees will be the sets of paths to the leaf nodes.

### Examples
    blocks = get_groups(tree)
"""
function get_groups(tree::AbstractTree; combine=0)
    leafnodes=get_leafnodes(tree)

    groups=Array{AbstractTree,1}[]
    nodes=AbstractTree[]

    for leaf in leafnodes
        node=history(leaf)[combine+1]
        if node ∉ nodes
            push!(nodes,node)
        end
    end

    for node in nodes
        scenario=history(node)
        children = collect(node)
        popfirst!(children)
        scenario = [scenario;children]
        push!(groups,scenario)
    end

    return groups
end

function get_scenarios(tree::AbstractTree)
    scenarios=get_groups(tree,combine=0)

    trees=AbstractTree[]

    for i in 1:length(scenarios)
        t=Leaf()
        t.ID=scenarios[i][1]
        t.name=t.ID.name
        for j in 2:length(scenarios[i])
            t=Tree([t])
            t.ID=scenarios[i][j]
            t.name=t.ID.name
        end
        push!(trees,t)
    end

    return trees
end

mutable struct Node
    children::Array{Node,1}
    pr::Float64
    name::String
end

function save_tree_to_file(tree::AbstractTree,filename::String)
    nodes=collect(tree,order=:breadth)

    f = open(filename,"w")
    println(f,"n,p")
    for n in nodes
        if parent(n)==nothing
            println(f,n.name*","*"-")
        else
            println(f,n.name*","*n.parent.name)
        end
    end
    close(f)
end

"""
	tree_from_leaves(leafnodes::Array{Array{Int64,1},1}, probs::Array{Float64,1})

Construct tree from Array of leaf nodes, and (optionally) the corresponding probabilities

### Required Arguments
`leafnodes` is an array of arrays defining the set of leaf nodes

### Optional Arguments
`probs` is an array of probabilities for the leaf nodes

### Example
    (tree,prob) = tree_from_leaves([[1,1,1],[1,1,2],[1,2,1],[1,2,2]],[0.25,0.25,0.25,0.25])
    tree = tree_from_leaves([[1,1,1],[1,1,2],[1,2,1],[1,2,2]])
"""#
function tree_from_leaves(leafnodes::Array{Array{Int64,1},1}, probs::Array{Float64,1})
    prob = Dict{AbstractTree,Float64}()

    function groupnode(node::Node)
        if length(node.children) == 0
            output = Leaf()
        else
            v = Array{AbstractTree,1}()
            for i = 1:length(node.children)
                push!(v, groupnode(node.children[i]))
                node.pr += node.children[i].pr
            end
            output = Tree(v)
            for c in output.children
                c.parent=output
            end
        end
        prob[output] = node.pr
        output
    end

    root = Node(Array{Node,1}(), 0.0, "")
    for i = 1:length(leafnodes)
        if leafnodes[i][1] != 1
            error("there must be a unique root node")
        end
        parent = root
        for j = 2:length(leafnodes[i])
            while leafnodes[i][j] > length(parent.children)
                n = Node(Array{Node,1}(), 0.0, "")
                push!(parent.children, n)
            end
            parent = parent.children[leafnodes[i][j]]
            if j == length(leafnodes[i])
                parent.pr = probs[i]
            end
        end
    end
    tree = groupnode(root)
    label_nodes(tree)
    return (tree, prob)
end

# Construct tree from Array of leaf nodes, without probabilities
function tree_from_leaves(leafnodes::Array{Array{Int64,1},1})
    function groupnode(node::Node)
        if length(node.children) == 0
            output = Leaf()
        else
            v = Array{AbstractTree,1}()
            for i = 1:length(node.children)
                push!(v, groupnode(node.children[i]))
            end
            output = Tree(v)
            for c in output.children
                c.parent=output
            end
        end
        output
    end

    root = Node(Array{Node,1}(), 0.0, "")
    for i = 1:length(leafnodes)
        if leafnodes[i][1] != 1
            error("there must be a unique root node")
        end
        parent = root
        for j = 2:length(leafnodes[i])
            while leafnodes[i][j] > length(parent.children)
                n = Node(Array{Node,1}(), 0.0, "")
                push!(parent.children, n)
            end
            parent = parent.children[leafnodes[i][j]]
            if j == length(leafnodes[i])
            end
        end
    end

    tree = groupnode(root)
    label_nodes(tree)
    return tree
end

# Construct tree from nested Vector, the first element of each vector is the
# probability of a node, remaining elements are vectors representing child nodes
function tree_from_nodes(nodes::Vector{Any})
    prob = Dict{AbstractTree,Float64}()

    function groupnode(node::Any, prob, prev)
        if length(node) == 1
            output = Leaf()
        else
            v = Array{AbstractTree,1}()
            for i = 2:length(node)
                push!(v, groupnode(node[i], prob, prev * node[1]))
            end
            output = Tree(v)
            for c in output.children
                c.parent=output
            end
        end
        prob[output] = node[1] * prev
        output
    end
    tree = groupnode(nodes, prob, 1.0)
    label_nodes(tree)
    return (tree, prob)
end

"""
	tree_from_file(filename::String)

Construct tree from a file, each line in the file is of the form B,A,...
representing an arc in the tree, from node "A" to node "B". The total
number of columns is arbitrary. The first row of the file should be
n,p,... these column headers are converted into symbols used to index
the `data`. Each column itself converted into a dictionary, indexed by
the node.

### Required Arguments
`string` is the full path of the file containing the tree

### Example
    tree, data = tree_from_file(joinpath(@__DIR__,"tree.csv"))
"""#
function tree_from_file(filename::String)
    data=Dict{Symbol,Dict{Node,Float64}}()
    data2=Dict{Symbol,Dict{AbstractTree,Float64}}()
    nodes=Dict{String,Node}()
    count=Dict{Node,Int64}()
    first=true
    f=open(filename)
    headers=Array{Symbol,1}()
    for l in eachline(f)
        if first
            first=false
            a=split(l,",")
            if a[1]!="n" || a[2]!="p"
                error("First row of tree file should be \"n,p,...\"")
            end
            for i in 3:length(a)
                push!(headers,Symbol(a[i]))
                data[Symbol(a[i])]=Dict{Node,Float64}()
                data2[Symbol(a[i])]=Dict{AbstractTree,Float64}()
            end
        else
            a=split(l,",")
            if length(a)==0
                @warn("Found blank line")
            elseif length(a)!=length(headers)+2
                @warn("Found line with wrong number of elements: \""*l*"\"")
                continue
            end

            if !((string)(a[2]) in keys(nodes))
                if (string)(a[2])!="-"
                    nodes[(string)(a[2])]=Node(Array{Node,1}(),1.0,(string)(a[2]))
                    count[nodes[(string)(a[2])]]=0
                end
            end
            if !((string)(a[1]) in keys(nodes))
                nodes[(string)(a[1])]=Node(Array{Node,1}(),1.0,(string)(a[1]))
                count[nodes[(string)(a[1])]]=0
            end
            if (string)(a[2])!="-"
                push!(nodes[(string)(a[2])].children,nodes[(string)(a[1])])
            end
            for i in 3:length(a)
                data[headers[i-2]][nodes[(string)(a[1])]]=parse(Float64, a[i])
            end
        end
    end

    close(f)

    for (nn,n) in nodes
        for i in n.children
            count[i]+=1
            if count[i]>1
                error("Node "*i.name*" defined multiple times")
            end
        end
    end

    found=0
    root=nothing
    for n in keys(count)
        if count[n]==0
            root=n
            found+=1
        end
    end

    if found==0
        error("No root node found")
    elseif found>1
        error("Multiple root nodes found")
    end

    function groupnode(node::Node)
        if length(node.children)==0
            output=Leaf(node.name)
        else
            v=Array{AbstractTree,1}()
            for n in node.children
                push!(v,groupnode(n))
            end
            output=Tree(v)
            output.name=node.name
            for c in output.children
                c.parent=output
            end
        end
        for h in headers
            data2[h][output]=data[h][node]
        end
        output
    end

    tree=groupnode(root)
    return (tree,data2)
end


text:
using Laszip

"""
Write XYZ data to a LIDAR laz (laszip compressed) or las format file. Usage:

	dat2las(FileName::AbstractString, xyz, hdr_vec=[]; scaleX=[], scaleY=[], scaleZ=[], offX=[], offY=[], offZ=[])

	Where:
		"FileName" Name of the output LIDAR file
		xyz  A Mx3 array with the point coordinates

Example. To write the x,y,z data to file "lixo.laz" do:

	dat2las("lixo.laz", xyz)
"""
function dat2las(fname::AbstractString, xyz, hdr_vec=[]; scaleX=[], scaleY=[], scaleZ=[], offX=[], offY=[], offZ=[])

	parse_inputs_dat2las(xyz, hdr_vec)
	n_rows, n_cols = size(xyz)

	if (!isempty(hdr_vec))
		min_x, max_x, min_y, max_y, min_z, max_z = hdr_vec[1:6]
	end

	#  Create the writer
	laszip_writer = convert(Ptr{Ptr{Cvoid}},pointer([pointer([0])]))
	if (laszip_create(laszip_writer) != 0)
		msgerror(laszip_writer, "creating laszip writer")
	end

	#header = pointer([pointer([create_empty_header()])])    # Get an empty header directly from C
	header = pointer([pointer([create_header()])])    # Get an empty header directly from C
	laszip_writer = unsafe_load(laszip_writer)

	if (laszip_get_header_pointer(laszip_writer, header) != 0)      # Get the header pointer
		msgerror(laszip_writer, "getting header pointer from laszip writer")
	end

	hdr = unsafe_load(unsafe_load(header))      # Get back the straight type

	# Populate the header
	hdr.version_major = 1
	hdr.version_minor = 2
	hdr.point_data_format = 1
	hdr.number_of_point_records = n_rows
	hdr.x_scale_factor = 1.0
	hdr.y_scale_factor = 1.0
	hdr.z_scale_factor = 1.0
	if (n_cols == 3 && isempty(hdr_vec))		# The 'regular' situation
		min_x, min_y, min_z = minimum(xyz, dims=1)
		max_x, max_y, max_z = maximum(xyz, dims=1)
	end
	hdr.max_x = max_x
	hdr.min_x = min_x
	hdr.max_y = max_y
	hdr.min_y = min_y
	hdr.max_z = max_z
	hdr.min_z = min_z

	# ----------------- Find reasonable scale_factor and offset -----------------------------------------
	if (hdr.min_x >= -360 && hdr.max_x <= 360 && hdr.min_y >= -90 && hdr.max_y <= 90)	# Assume geogs
		if (isempty(scaleX))  hdr.x_scale_factor = 1e-7
		else                  hdr.x_scale_factor = scaleX
		end
		if (isempty(scaleY))  hdr.y_scale_factor = 1e-7
		else                  hdr.y_scale_factor = scaleY
		end
	else
		if (isempty(scaleY))  hdr.x_scale_factor = 1e-3
		else                  hdr.x_scale_factor = scaleX
		end
		if (isempty(scaleY))  hdr.y_scale_factor = 1e-3
		else                  hdr.y_scale_factor = scaleY
		end
	end
	if (isempty(scaleZ))  hdr.z_scale_factor = 1e-2
	else                  hdr.z_scale_factor = scaleZ
	end

	if (!isnan(hdr.min_x) && !isnan(hdr.max_x))
		if (isempty(offX))
			hdr.x_offset = (floor((hdr.min_x + hdr.max_x)/hdr.x_scale_factor/20000000)) * 10000000 * hdr.x_scale_factor
		else
			hdr.x_offset = offX
		end
	end
	if (!isnan(hdr.min_y) && !isnan(hdr.max_y))
		if (isempty(offY))
			hdr.y_offset = (floor((hdr.min_y + hdr.max_y)/hdr.y_scale_factor/20000000)) * 10000000 * hdr.y_scale_factor
		else
			hdr.y_offset = offY
		end
	end
	if (!isnan(hdr.min_z) && !isnan(hdr.max_z))
		if (isempty(offZ))
			hdr.z_offset = (floor((hdr.min_z + hdr.max_z)/hdr.z_scale_factor/20000000)) * 10000000 * hdr.z_scale_factor
		else
			hdr.z_offset = offZ
		end
	end
	# ---------------------------------------------------------------------------------------------------

	# This is the case where we are storing a grid pretending it's a regular LAZ file. Must hijack some header members
	if (!isempty(hdr_vec))
		hdr.x_scale_factor = hdr.z_scale_factor		# Because in fact we only have zz's
		hdr.y_scale_factor = hdr.z_scale_factor
		hdr.x_offset = hdr.z_offset
		hdr.y_offset = hdr.z_offset
		hdr.number_of_point_records = UInt32(ceil(n_rows / 3))
		hdr.global_encoding = 32768;	# Use this number to codify as GRID. bin(UInt16(32768)) = "1000000000000000"

		one = (hdr_vec[7] == 0 ? 1 : 0)
		hdr.project_ID_GUID_data_1 = hdr_vec[7]
		hdr.project_ID_GUID_data_2 = round(Int, (hdr.max_y - hdr.min_y) / hdr_vec[8]) + one 	# n_rows in 2D array
		hdr.project_ID_GUID_data_3 = round(Int, (hdr.max_x - hdr.min_x) / hdr_vec[9]) + one 	# n_cols in 2D array

		#hdr.system_identifier.d1 = UInt8('G')	# This is ugly. There must be a better way
		#hdr.system_identifier.d2 = UInt8('R')
		#hdr.system_identifier.d3 = UInt8('D')
		#mutateit(hdr.system_identifier, "d1", 'G')
		#mutateit(hdr.system_identifier, "d2", 'R')
		#mutateit(hdr.system_identifier, "d3", 'D')
	end

	# Save back the header to its C pointer
	unsafe_store!(unsafe_load(header), hdr)

	if (laszip_open_writer(laszip_writer, fname, 1) != 0)
		msgerror(laszip_writer, @sprintf("opening laszip writer for %s", fname))
	end

	# Get a pointer to the points that will be written
	#point = pointer([pointer([create_empty_point()])])
	point = pointer([pointer([create_point()])])
	if (laszip_get_point_pointer(laszip_writer, point) != 0)
		msgerror(laszip_writer, "getting point pointer from laszip writer")
	end
	point = unsafe_load(point)

	coordinates = zeros(3)
	if (n_cols == 3)
		for n = 1:n_rows
			#=
			pt = unsafe_load(point)
			pt.X = round(Int32, (xyz[n,1]-hdr.x_offset) * hdr.x_scale_factor)
			pt.Y = round(Int32, (xyz[n,2]-hdr.y_offset) * hdr.y_scale_factor)
			pt.Z = round(Int32, (xyz[n,3]-hdr.z_offset) * hdr.z_scale_factor)
			unsafe_store!(point, pt)
			=#
			coordinates[1] = xyz[n,1]
			coordinates[2] = xyz[n,2]
			coordinates[3] = xyz[n,3]
			laszip_set_coordinates(laszip_writer, convert(Ptr{Cdouble}, pointer(coordinates)))
			laszip_write_point(laszip_writer)
		end
	else
		# Deal with special case of a Z column only where we will cheat by splitting trhough the XYZ
		r = rem(n_rows, 3)
		last_ind = (r == 0 ? n_rows : n_rows - 1)
		for n = 1:3:last_ind
			#pt = unsafe_load(point)
			#pt.X = round(Int32, (xyz[n,1]-hdr.z_offset)   * hdr.z_scale_factor)
			#pt.Y = round(Int32, (xyz[n+1,1]-hdr.z_offset) * hdr.z_scale_factor)
			#pt.Z = round(Int32, (xyz[n+2,1]-hdr.z_offset) * hdr.z_scale_factor)
			coordinates[1] = xyz[n]
			coordinates[2] = xyz[n+1]
			coordinates[3] = xyz[n+2]
			laszip_set_coordinates(laszip_writer, convert(Ptr{Cdouble}, pointer(coordinates)))
			laszip_write_point(laszip_writer)
		end
		#pt = unsafe_load(point)
		if (r == 1)
			#pt.X = round(Int32, (xyz[last_ind+1,1]-hdr.z_offset) * hdr.z_scale_factor)
			#pt.Y = pt.X
			#pt.Z = pt.X
			coordinates[1] = xyz[last_ind+1]
			laszip_set_coordinates(laszip_writer, convert(Ptr{Cdouble}, pointer(coordinates)))
			laszip_write_point(laszip_writer)
		elseif (r == 2)
			#pt.X = round(Int32, (xyz[last_ind+1,1]-hdr.z_offset) * hdr.z_scale_factor)
			#pt.Y = round(Int32, (xyz[last_ind+2,1]-hdr.z_offset) * hdr.z_scale_factor)
			#pt.Z = pt.Y
			coordinates[1] = xyz[last_ind+1]
			coordinates[2] = xyz[last_ind+2]
			laszip_set_coordinates(laszip_writer, convert(Ptr{Cdouble}, pointer(coordinates)))
			laszip_write_point(laszip_writer)
		end
	end

	# Close the writer
	if (laszip_close_writer(laszip_writer) != 0)
		msgerror(laszip_writer, "closing laszip writer")
	end
	# Destroy the writer
	if (laszip_destroy(laszip_writer) != 0)
		msgerror(laszip_writer, "destroying laszip writer")
	end

end

# --------------------------------------------------------------------------
function parse_inputs_dat2las(xyz, hdr_vec=[])
# Check validity of input and in future will parse string options

	n_rows, n_cols = size(xyz)
	if ((n_cols != 3 && n_cols != 1))
		error("Input array can only have 1 or 3 columns OR be a 2D array")
	end
	if (!isempty(hdr_vec))
		if (length(hdr_vec) < 7)
			error("HDR argument does not have at least 7 elements")
		end
	end
end

ext: jl

max_stars_repo_name: joa-quim/Laszip

max_issues_repo_name: joa-quim/Laszip

max_forks_repo_name: joa-quim/Laszip


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Runge-Kutta 4th order solver
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

struct SimpleRK4 <: DiffEqBase.AbstractODEAlgorithm end
export SimpleRK4

mutable struct SimpleRK4Integrator{IIP, S, T, P, F} <: DiffEqBase.AbstractODEIntegrator{SimpleRK4, IIP, S, T}
    f::F             # ..................................... Equations of motion
    uprev::S         # .......................................... Previous state
    u::S             # ........................................... Current state
    tmp::S           #  Auxiliary variable similar to state to avoid allocations
    tprev::T         # ...................................... Previous time step
    t::T             # ....................................... Current time step
    t0::T            # ........... Initial time step, only for re-initialization
    dt::T            # ............................................... Step size
    tdir::T          # ........................................ Not used for RK4
    p::P             # .................................... Parameters container
    u_modified::Bool # ..... If `true`, then the input of last step was modified
    ks::Vector{S}    # ........................... Interpolants of the algorithm
end

const SRK4 = SimpleRK4Integrator

# If `true`, then the equation of motion format is `f!(du,u,p,t)` instead of
# `du = f(u,p,t)`.
DiffEqBase.isinplace(::SRK4{IIP}) where {IIP} = IIP

################################################################################
#                                Initialization
################################################################################

function DiffEqBase.__init(prob::ODEProblem, alg::SimpleRK4;
                           dt = error("dt is required for this algorithm"))
    simplerk4_init(prob.f,
                   DiffEqBase.isinplace(prob),
                   prob.u0,
                   prob.tspan[1],
                   dt,
                   prob.p)
end

function DiffEqBase.__solve(prob::ODEProblem, alg::SimpleRK4;
                            dt = error("dt is required for this algorithm"))
    u0    = prob.u0
    tspan = prob.tspan
    ts    = Array(tspan[1]:dt:tspan[2])
    n     = length(ts)
    us    = Vector{typeof(u0)}(undef, n)

    @inbounds us[1] = _copy(u0)

    integ = simplerk4_init(prob.f, DiffEqBase.isinplace(prob), prob.u0,
                           prob.tspan[1], dt, prob.p)

    # FSAL
    for i = 1:n-1
        step!(integ)
        us[i+1] = _copy(integ.u)
    end

    sol = DiffEqBase.build_solution(prob, alg, ts, us, calculate_error = false)

    DiffEqBase.has_analytic(prob.f) &&
        DiffEqBase.calculate_solution_errors!(sol;
                                              timeseries_errors = true,
                                              dense_errors = false)

    return sol
end

@inline function simplerk4_init(f::F, IIP::Bool, u0::S, t0::T, dt::T, p::P) where
    {F, P, T, S<:AbstractArray{T}}

    # Allocate the vector with the interpolants. For RK4, we need 5.
    ks = [zero(u0) for _ = 1:5]

    integ = SRK4{IIP, S, T, P, F}(f,
                                  _copy(u0),
                                  _copy(u0),
                                  _copy(u0),
                                  t0,
                                  t0,
                                  t0,
                                  dt,
                                  sign(dt),
                                  p,
                                  true,
                                  ks)

    return integ
end

################################################################################
#                                   Stepping
################################################################################

@inline @muladd function DiffEqBase.step!(integ::SRK4{true, S, T}) where {T, S}
    integ.uprev       .= integ.u
    tmp                = integ.tmp
    f!                 = integ.f
    p                  = integ.p
    t                  = integ.t
    dt                 = integ.dt
    uprev              = integ.uprev
    u                  = integ.u
    k₁, k₂, k₃, k₄, k₅ = integ.ks

    # If the input was modified, then we need to recompute the initial pass of
    # this step.
    if integ.u_modified
        f!(k₁, uprev, p, t)
        integ.u_modified = false
    else
        k₁ .= k₅
    end

    # This algorithm is faster than the one using broadcasts.
    @inbounds begin
        L = length(u)

        for i in 1:L
            tmp[i] = uprev[i] + dt*k₁[i]/2
        end
        f!(k₂, tmp, p, t + dt/2)

        for i in 1:L
            tmp[i] = uprev[i] + dt*k₂[i]/2
        end
        f!(k₃, tmp, p, t + dt/2)

        for i in 1:L
            tmp[i] = uprev[i] + dt*k₃[i]
        end
        f!(k₄, tmp, p, t + dt)

        for i = 1:L
            u[i] = uprev[i] + (dt/6)*( 2*(k₂[i] + k₃[i]) + (k₁[i] + k₄[i]) )
        end

        f!(k₅, u, p, t + dt)
    end

    integ.tprev = t
    integ.t += dt

    return nothing
end

@inline @muladd function DiffEqBase.step!(integ::SRK4{false, S, T}) where {T, S}
    integ.uprev = integ.u
    f           = integ.f
    p           = integ.p
    t           = integ.t
    dt          = integ.dt
    uprev       = integ.uprev

    # If the input was modified, then we need to recompute the initial pass of
    # this step.
    if integ.u_modified
        k₁ = f(uprev, p, t)
        integ.u_modified = false
    else
        @inbounds k₁ = integ.ks[5]
    end

    tmp = uprev + dt*k₁/2
    k₂ = f(tmp, p, t + dt/2)

    tmp = uprev + dt*k₂/2
    k₃ = f(tmp, p, t + dt/2)

    tmp = uprev + dt*k₃
    k₄ = f(tmp, p, t + dt)

    integ.u = uprev + (dt/6)*( 2*(k₂ + k₃) + (k₁ + k₄) )
    k₅ = f(integ.u, p, t + dt)

    # Update the interpolants in the integrator. This is necessary for the
    # interpolation.
    @inbounds begin
        integ.ks[1] = k₁
        integ.ks[2] = k₂
        integ.ks[3] = k₃
        integ.ks[4] = k₄
        integ.ks[5] = k₅
    end

    integ.tprev = t
    integ.t += dt

    return nothing
end

################################################################################
#                                Interpolation
################################################################################

@inline @muladd function (integ::SRK4)(t::T) where T
    t₁, t₀, dt = integ.t, integ.tprev, integ.dt

    y₀ = integ.uprev
    y₁ = integ.u
    ks = integ.ks
    Θ  = (t - t₀)/dt

    # Hermite interpolation.
    @inbounds if !isinplace(integ)
        u = (1-Θ)*y₀ + Θ*y₁ + Θ*(Θ-1)*( (1-2Θ)*(y₁-y₀) +
                                        (Θ-1)*dt*ks[1] +
                                        Θ*dt*ks[5])
        return u
    else
        u = similar(y₁)
        for i in 1:length(u)
            u[i] = (1-Θ)*y₀[i] + Θ*y₁[i] + Θ*(Θ-1)*( (1-2Θ)*(y₁[i]-y₀[i])+
                                                     (Θ-1)*dt*ks[1][i] +
                                                     Θ*dt*ks[5][i])
        end

        return u
    end
end


text:
#=
# Import a simple tet or hex mesh from a .msh file.
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("unstruct3d");

useLog("unstruct3dlog", level=3)

domain(3, grid=UNSTRUCTURED)
functionSpace(order=2)

# Unit cube
mesh("utet.msh")  # Using tets
#mesh("uhex.msh")   # Using hexes

u = variable("u")
testSymbol("v")

boundary(u, 1, DIRICHLET, "sin(0.5*pi*x)*sin(pi*y)*sin(pi*z)")

coefficient("K", "x*x + 4")
coefficient("C", 10)
coefficient("f", "x*pi*cos(pi*0.5*x)*sin(pi*y)*sin(pi*z) + (-10-(x*x+4)*pi*pi*2.25)*sin(pi*0.5*x)*sin(pi*y)*sin(pi*z)")

# PDE: div(K*grad(u)) -C*u = f
weakForm(u, "K*dot(grad(u), grad(v)) + C*u*v + f*v")

solve(u);

# exact solution is sin(pi*0.5*x)*sin(pi*y)*sin(pi*z)
# check error
maxerr = 0;
allerr = zeros(size(u.values));
exact(x,y,z) = sin(0.5*pi*x)*sin(pi*y)*sin(pi*z);

for i=1:size(Finch.grid_data.allnodes,2)
    x = Finch.grid_data.allnodes[1,i];
    y = Finch.grid_data.allnodes[2,i];
    z = Finch.grid_data.allnodes[3,i];
    err = abs(u.values[i] - exact(x,y,z));
    allerr[i] = err;
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

# using Plots
# pyplot();
# display(plot(Finch.grid_data.allnodes[1,:], Finch.grid_data.allnodes[2,:], allerr[:], st=:surface))

finalize_finch();

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
module test_Plot7

import ..ModiaMath
using Test
t = range(0.0, stop=10.0, length=100)

# Define signals
mutable struct MyThermodynamicState
    p::Float64
    T::Float64
end

specificEnthalpy(state::MyThermodynamicState) = 2.0*state.T
dynamicViscosity(state::MyThermodynamicState) = 2.0*state.p

const dependentVariables = Dict{AbstractString,Function}("h"   => specificEnthalpy,
                                                         "eta" => dynamicViscosity)

ModiaMath.variablesDependingOnStruct(state::MyThermodynamicState) = dependentVariables

state = MyThermodynamicState[]
for i = 1:length(t)
    push!(state, MyThermodynamicState(i*2.0,i*3.0))
end

result = Dict{AbstractString,Any}("time" => t, "state" => state)

# Plots
ModiaMath.plot(result, ("state.p", "state.T", "state.h", "state.eta"), figure=7)

end


text:
# generic APIs

## dense vector descriptor

mutable struct CuDenseVectorDescriptor
    handle::cusparseDnVecDescr_t

    function CuDenseVectorDescriptor(x::DenseCuVector)
        desc_ref = Ref{cusparseDnVecDescr_t}()
        cusparseCreateDnVec(desc_ref, length(x), x, eltype(x))
        obj = new(desc_ref[])
        finalizer(cusparseDestroyDnVec, obj)
        obj
    end
end

Base.unsafe_convert(::Type{cusparseDnVecDescr_t}, desc::CuDenseVectorDescriptor) = desc.handle

## sparse vector descriptor

mutable struct CuSparseVectorDescriptor
    handle::cusparseSpVecDescr_t

    function CuSparseVectorDescriptor(x::CuSparseVector, IndexBase::Char)
        desc_ref = Ref{cusparseSpVecDescr_t}()
        cusparseCreateSpVec(desc_ref, length(x), nnz(x), nonzeroinds(x), nonzeros(x),
                            eltype(nonzeroinds(x)), IndexBase, eltype(x))
        obj = new(desc_ref[])
        finalizer(cusparseDestroySpVec, obj)
        obj
    end
end

Base.unsafe_convert(::Type{cusparseSpVecDescr_t}, desc::CuSparseVectorDescriptor) = desc.handle

## dense matrix descriptor

mutable struct CuDenseMatrixDescriptor
    handle::cusparseDnMatDescr_t

    function CuDenseMatrixDescriptor(x::DenseCuMatrix)
        desc_ref = Ref{cusparseDnMatDescr_t}()
        cusparseCreateDnMat(desc_ref, size(x)..., stride(x,2), x, eltype(x), CUSPARSE_ORDER_COL)
        obj = new(desc_ref[])
        finalizer(cusparseDestroyDnMat, obj)
        obj
    end
end

Base.unsafe_convert(::Type{cusparseDnMatDescr_t}, desc::CuDenseMatrixDescriptor) = desc.handle

## sparse matrix descriptor

mutable struct CuSparseMatrixDescriptor
    handle::cusparseSpMatDescr_t

    function CuSparseMatrixDescriptor(A::CuSparseMatrixCSR, IndexBase::Char)
        desc_ref = Ref{cusparseSpMatDescr_t}()
        cusparseCreateCsr(
            desc_ref,
            size(A)..., length(nonzeros(A)),
            A.rowPtr, A.colVal, nonzeros(A),
            eltype(A.rowPtr), eltype(A.colVal), IndexBase, eltype(nonzeros(A))
        )
        obj = new(desc_ref[])
        finalizer(cusparseDestroySpMat, obj)
        return obj
    end

    function CuSparseMatrixDescriptor(A::CuSparseMatrixCSC, IndexBase::Char; convert=true)
        desc_ref = Ref{cusparseSpMatDescr_t}()
        if convert
            # many algorithms, e.g. mv! and mm!, do not support CSC sparse format
            # so we eagerly convert this to a CSR matrix.
            cusparseCreateCsr(
                desc_ref,
                reverse(size(A))..., length(nonzeros(A)),
                A.colPtr, rowvals(A), nonzeros(A),
                eltype(A.colPtr), eltype(rowvals(A)), IndexBase, eltype(nonzeros(A))
            )
        else
            cusparseCreateCsc(
                desc_ref,
                size(A)..., length(nonzeros(A)),
                A.colPtr, rowvals(A), nonzeros(A),
                eltype(A.colPtr), eltype(rowvals(A)), IndexBase, eltype(nonzeros(A))
            )
        end
        obj = new(desc_ref[])
        finalizer(cusparseDestroySpMat, obj)
        return obj
    end
end

Base.unsafe_convert(::Type{cusparseSpMatDescr_t}, desc::CuSparseMatrixDescriptor) = desc.handle

## API functions

function gather!(X::CuSparseVector, Y::CuVector, index::SparseChar)
    descX = CuSparseVectorDescriptor(X, index)
    descY = CuDenseVectorDescriptor(Y)

    cusparseGather(handle(), descY, descX)

    X
end

function mv!(transa::SparseChar, alpha::Number, A::Union{CuSparseMatrixBSR{T},CuSparseMatrixCSR{T}},
             X::DenseCuVector{T}, beta::Number, Y::DenseCuVector{T}, index::SparseChar) where {T}
    m,n = size(A)

    if transa == 'N'
        chkmvdims(X,n,Y,m)
    elseif transa == 'T' || transa == 'C'
        chkmvdims(X,m,Y,n)
    end

    descA = CuSparseMatrixDescriptor(A, index)
    descX = CuDenseVectorDescriptor(X)
    descY = CuDenseVectorDescriptor(Y)
    compute_type = T == Float16 && version() >= v"11.4.0" ? Float32 : T

    function bufferSize()
        out = Ref{Csize_t}()
        cusparseSpMV_bufferSize(handle(), transa, Ref{compute_type}(alpha), descA, descX, Ref{compute_type}(beta),
                                descY, compute_type, CUSPARSE_SPMV_ALG_DEFAULT, out)
        return out[]
    end
    with_workspace(bufferSize) do buffer
        cusparseSpMV(handle(), transa, Ref{compute_type}(alpha), descA, descX, Ref{compute_type}(beta),
                     descY, compute_type, CUSPARSE_SPMV_ALG_DEFAULT, buffer)
    end
    Y
end

function mv!(transa::SparseChar, alpha::Number, A::CuSparseMatrixCSC{T}, X::DenseCuVector{T},
             beta::Number, Y::DenseCuVector{T}, index::SparseChar) where {T}
    ctransa = 'N'
    if transa == 'N'
        ctransa = 'T'
    elseif transa == 'C' && T <: Complex
        throw(ArgumentError("Matrix-vector multiplication with the adjoint of a CSC matrix" *
                            " is not supported. Use a CSR matrix instead."))
    end

    n,m = size(A)

    if ctransa == 'N'
        chkmvdims(X,n,Y,m)
    elseif ctransa == 'T' || ctransa == 'C'
        chkmvdims(X,m,Y,n)
    end

    descA = CuSparseMatrixDescriptor(A, index)
    descX = CuDenseVectorDescriptor(X)
    descY = CuDenseVectorDescriptor(Y)
    compute_type = T == Float16 && version() >= v"11.4.0" ? Float32 : T

    function bufferSize()
        out = Ref{Csize_t}()
        cusparseSpMV_bufferSize(handle(), ctransa, Ref{compute_type}(alpha), descA, descX, Ref{compute_type}(beta),
                                descY, compute_type, CUSPARSE_SPMV_ALG_DEFAULT, out)
        return out[]
    end
    with_workspace(bufferSize) do buffer
        cusparseSpMV(handle(), ctransa, Ref{compute_type}(alpha), descA, descX, Ref{compute_type}(beta),
                     descY, compute_type, CUSPARSE_SPMV_ALG_DEFAULT, buffer)
    end

    return Y
end

function mm!(transa::SparseChar, transb::SparseChar, alpha::Number, A::CuSparseMatrixCSR{T},
             B::DenseCuMatrix{T}, beta::Number, C::DenseCuMatrix{T}, index::SparseChar) where {T}
    m,k = size(A)
    n = size(C)[2]

    if transa == 'N' && transb == 'N'
        chkmmdims(B,C,k,n,m,n)
    elseif transa == 'N' && transb != 'N'
        chkmmdims(B,C,n,k,m,n)
    elseif transa != 'N' && transb == 'N'
        chkmmdims(B,C,m,n,k,n)
    elseif transa != 'N' && transb != 'N'
        chkmmdims(B,C,n,m,k,n)
    end

    descA = CuSparseMatrixDescriptor(A, index)
    descB = CuDenseMatrixDescriptor(B)
    descC = CuDenseMatrixDescriptor(C)

    function bufferSize()
        out = Ref{Csize_t}()
        cusparseSpMM_bufferSize(
            handle(), transa, transb, Ref{T}(alpha), descA, descB, Ref{T}(beta),
            descC, T, CUSPARSE_SPMM_ALG_DEFAULT, out)
        return out[]
    end
    with_workspace(bufferSize) do buffer
        cusparseSpMM(
            handle(), transa, transb, Ref{T}(alpha), descA, descB, Ref{T}(beta),
            descC, T, CUSPARSE_SPMM_ALG_DEFAULT, buffer)
    end

    return C
end

function mm!(transa::SparseChar, transb::SparseChar, alpha::Number, A::CuSparseMatrixCSC{T},
             B::DenseCuMatrix{T}, beta::Number, C::DenseCuMatrix{T}, index::SparseChar) where {T}
    ctransa = 'N'
    if transa == 'N'
        ctransa = 'T'
    elseif transa == 'C' && T <: Complex
        throw(ArgumentError("Matrix-matrix multiplication with the adjoint of a CSC matrix" *
                            " is not supported. Use a CSR matrix instead."))
    end

    k,m = size(A)
    n = size(C)[2]

    if ctransa == 'N' && transb == 'N'
        chkmmdims(B,C,k,n,m,n)
    elseif ctransa == 'N' && transb != 'N'
        chkmmdims(B,C,n,k,m,n)
    elseif ctransa != 'N' && transb == 'N'
        chkmmdims(B,C,m,n,k,n)
    elseif ctransa != 'N' && transb != 'N'
        chkmmdims(B,C,n,m,k,n)
    end

    descA = CuSparseMatrixDescriptor(A, index)
    descB = CuDenseMatrixDescriptor(B)
    descC = CuDenseMatrixDescriptor(C)

    function bufferSize()
        out = Ref{Csize_t}()
        cusparseSpMM_bufferSize(
            handle(), ctransa, transb, Ref{T}(alpha), descA, descB, Ref{T}(beta),
            descC, T, CUSPARSE_SPMM_ALG_DEFAULT, out)
        return out[]
    end
    with_workspace(bufferSize) do buffer
        cusparseSpMM(
            handle(), ctransa, transb, Ref{T}(alpha), descA, descB, Ref{T}(beta),
            descC, T, CUSPARSE_SPMM_ALG_DEFAULT, buffer)
    end

    return C
end

function mm!(transa::SparseChar, transb::SparseChar, α::Number, A::CuSparseMatrixCSR{T},
             B::CuSparseMatrixCSR{T}, β::Number, C::CuSparseMatrixCSR{T}, index::SparseChar) where {T}
    m,k = size(A)
    n = size(C)[2]
    alpha = convert(T, α)
    beta  = convert(T, β)

    if transa == 'N' && transb == 'N'
        chkmmdims(B,C,k,n,m,n)
    else
        throw(ArgumentError("Sparse mm! only supports transa ($transa) = 'N' and transb ($transb) = 'N'"))
    end

    descA = CuSparseMatrixDescriptor(A, index)
    descB = CuSparseMatrixDescriptor(B, index)
    descC = CuSparseMatrixDescriptor(C, index)

    spgemm_Desc = CuSpGEMMDescriptor()
    function buffer1Size()
        out = Ref{Csize_t}(0)
        cusparseSpGEMM_workEstimation(
            handle(), transa, transb, Ref{T}(alpha), descA, descB, Ref{T}(beta),
            descC, T, CUSPARSE_SPGEMM_DEFAULT, spgemm_Desc, out, CU_NULL)
        return out[]
    end
    with_workspace(buffer1Size; keep=true) do buffer
        out = Ref{Csize_t}(sizeof(buffer))
        cusparseSpGEMM_workEstimation(
            handle(), transa, transb, Ref{T}(alpha), descA, descB, Ref{T}(beta),
            descC, T, CUSPARSE_SPGEMM_DEFAULT, spgemm_Desc, out, buffer)
    end
    function buffer2Size()
        out = Ref{Csize_t}(0)
        cusparseSpGEMM_compute(
            handle(), transa, transb, Ref{T}(alpha), descA, descB, Ref{T}(beta),
            descC, T, CUSPARSE_SPGEMM_DEFAULT, spgemm_Desc, out, CU_NULL)
        return out[]
    end
    with_workspace(buffer2Size; keep=true) do buffer
        out = Ref{Csize_t}(sizeof(buffer))
        cusparseSpGEMM_compute(
            handle(), transa, transb, Ref{T}(alpha), descA, descB, Ref{T}(beta),
            descC, T, CUSPARSE_SPGEMM_DEFAULT, spgemm_Desc, out, buffer)
    end
    Cm   = Ref{Int64}()
    Cn   = Ref{Int64}()
    Cnnz1 = Ref{Int64}()
    cusparseSpMatGetSize(descC, Cm, Cn, Cnnz1)
    # SpGEMM_copy assumes A*B and C have the same sparsity pattern if
    # beta is not zero. If that isn't the case, we must use broadcasted
    # add to get the correct result.
    if beta == zero(beta)
        unsafe_free!(C.rowPtr)
        unsafe_free!(C.colVal)
        unsafe_free!(C.nzVal)
        C.rowPtr = CuVector{Cint}(undef, Cm[] + 1)
        C.colVal = CuVector{Cint}(undef, Cnnz1[])
        C.nzVal  = CuVector{T}(undef, Cnnz1[])
        C.nnz    = Cnnz1[]
        cusparseCsrSetPointers(descC, C.rowPtr, C.colVal, C.nzVal)
        cusparseSpGEMM_copy(handle(), transa, transb, Ref{T}(alpha), descA, descB,
                            Ref{T}(beta), descC, T, CUSPARSE_SPGEMM_DEFAULT, spgemm_Desc)
        return C
    else
        newC = CuSparseMatrixCSR(CUDA.zeros(Cint, Cm[] + 1), CUDA.zeros(Cint, Cnnz1[]), CUDA.zeros(T, Cnnz1[]), (Cm[], Cn[]))
        descnewC = CuSparseMatrixDescriptor(newC, index)
        cusparseSpGEMM_copy(handle(), transa, transb, Ref{T}(alpha), descA, descB,
                            Ref{T}(beta), descnewC, T, CUSPARSE_SPGEMM_DEFAULT, spgemm_Desc)
        D = beta.*C .+ newC
        return D
    end
end


text:
using Test

using QuantumCircuits
using QuantumCircuits.QCircuits.Math

using QuantumCircuits.Simulation.Gates

gates = [ZZ, XX, YY]
θs = [0, π/4, π/2, 3π/4, π, 5π/4, 3π/2, 7π/4, 2π]

for g in gates
    for θ in θs
        simqc = QCircuit(2)
        g(simqc, 0, 1, θ, true)
        m1 = tomatrix(simqc)

        simqc = QCircuit(2)
        g(simqc, 0, 1, θ, false)
        m2 = tomatrix(simqc)

        @test unitary_error(m1, m2) < 1e-8
    end
end


text:
### A Pluto.jl notebook ###
# v0.18.2

using Markdown
using InteractiveUtils

# ╔═╡ e506c710-a422-11ec-08e1-cd6c187d3881
using CommonMark

# ╔═╡ 8d5d5caa-1012-4e4e-af93-9373bf8eac71

# ╔═╡ 362039ef-eee1-4c1a-b1f9-b95b803d816b
cm"""
# Simple shit not working
"""

# ╔═╡ 8b295661-4a95-451f-97ae-1c8100de8adc
cm"""
```math
x = \frac{\sqrt{x^2-1}}{32}
```
"""

# ╔═╡ 582ac6a4-9e64-4021-b322-157c9ba7013f
cm"""
```math
\ce{H2SO4}
```
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CommonMark = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"

[compat]
CommonMark = "~0.8.6"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "4cd7063c9bdebdbd55ede1af70f3c2f48fab4215"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.6"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
"""

# ╔═╡ Cell order:
# ╠═e506c710-a422-11ec-08e1-cd6c187d3881
# ╠═8d5d5caa-1012-4e4e-af93-9373bf8eac71
# ╟─362039ef-eee1-4c1a-b1f9-b95b803d816b
# ╟─8b295661-4a95-451f-97ae-1c8100de8adc
# ╟─582ac6a4-9e64-4021-b322-157c9ba7013f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: mdAshford/youCANdoThermodynamics

max_issues_repo_name: mdAshford/youCANdoThermodynamics

max_forks_repo_name: mdAshford/youCANdoThermodynamics


text:
using ExprParsers

# ExprParsers.to_expr(x::Union{DataType, UnionAll}) = Base.Meta.parse(repr(x))
function ExprParsers.to_expr(tv::TypeVar)
  if tv.lb === Union{} && tv.ub === Any
    tv.name
  elseif tv.lb === Union{}
    :($(tv.name) <: $(tv.ub))
  elseif tv.ub === Any
    :($(tv.name) >: $(tv.lb))
  else
    :($(tv.lb) <: $(tv.name) <: $(tv.ub))
  end
end


text:
using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
using PGFPlotsX, BSON, CodecZstd, DataFrames

push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"""
\usepackage{xcolor}
\definecolor{mblue}{HTML}{0992F2}
\definecolor{morange}{HTML}{ff7f0e}
\definecolor{mgreen}{HTML}{2ca02c}
\definecolor{mred}{HTML}{d62728}
""")

include(joinpath(@__DIR__, "paths.jl"))

results = open(joinpath(datapath, "results.bson.zstd")) do f
    s = ZstdDecompressorStream(f)
    res = BSON.load(s)[:results]
    close(s)
    res
end

include(joinpath(@__DIR__, "order.jl"))

sort!(results, :logp_hat)
sort!(results, :experiment, by = x -> findfirst(==(x), experiment_order))
sort!(results, :model, by = x -> findfirst(==(x), model_order))
best = combine(groupby(results, [:model, :experiment]),
               names(results) .=> x -> [last(x)], renamecols = false)
best.logp_hat_rel = best.logp_hat .- repeat(best.logp_hat[best.model .== :Baseline],
                                            length(union(best.model)))

groups = groupby(best[best.model .!= :Baseline, :], :model)
y = [g.logp_hat_rel for g in groups]

function toname(e)
    e == :Raby07_planningforbreakfast && return "Raby07 planning"
    e == :Raby07_breakfastchoice && return "Raby07 choice"
    replace(string(e), "_" => " ")
end

f1 = @pgf Axis({legend_entries =["Replay-and-Plan", "Plastic Caching",
                                 "Episodic-Like Memory", "Motivational Control"],
           legend_columns = 2,
           font = raw"\scriptsize", axis_lines = "left",
           x_axis_line_style = "white", ymax = 150, ymin = -2,
           y_axis_line_style = "-",
           "every axis legend/.append style" = {at = "{(0.01, 1)}",
                                                anchor = "south west"},
           xtick = 1:28, xmin = .5, xmax = 23.5,
           cycle_list = ["{thick, morange, mark = square}",
                         "{thick, mred, mark = square}",
                         "{thick, mblue, mark = square}",
                         "{thick, mgreen, mark = square}"],
           xticklabels = toname.(groups[1].experiment),
           width = "9cm", height = "6cm",
           ylabel = raw"$\Delta \log\hat p$",
           xticklabel_style = {rotate = 60, anchor = "east"}},
           "\\draw[draw = none, fill = mred!12] (-.5, 0) rectangle (11.5, 150);",
           "\\draw[draw = none, fill = mblue!12] (11.5, 0) rectangle (17.5, 150);",
           "\\draw[draw = none, fill = mgreen!12] (17.5, 0) rectangle (23, 150);",
           raw"\node[anchor = center] at (6, 135) {\color{black}planning experiments};",
           raw"\node[anchor = center] at (14.5, 135) {\parbox{2cm}{\centering memory experiments}};",
           raw"\node[anchor = center] at (20, 135) {\parbox{2cm}{\centering satiety experiments}};",
          [Plot(Coordinates(1:length(yi), yi)) for yi in y]...,
           Plot({no_marks, dashed, domain = "{0:28}"}, Expression("0")),
         )
pgfsave(joinpath(figpath, "logp_hat_rel_experiments.tikz"), f1)

pe = 1:11
me = 12:17
se = 18:22
ce = 23
ys = [sum(yi[e]) for yi in y, e in (ce, se, me, pe)]
f2 = @pgf Axis({ybar_stacked, xtick = 1:4, font = raw"\scriptsize", axis_lines =
                "left", x_axis_line_style = "white", xmin = 0.5, xmax = 4.5,
                height = "6cm", width = "6cm", ymax = 820, ymin = 0,
                y_axis_line_style = "-", "thin",
                cycle_list = ["fill = gray", "fill = mgreen!12",
                              "fill = mblue!12", "fill = mred!12"],
                ylabel = raw"$\Delta \log\hat p$",
                xticklabels =
                [raw"\parbox{3cm}{\color{morange}\raggedleft Replay-and-Plan\\\&
                                  Episodic-Like Memory\\\& Motivational Control}",
                 raw"\parbox{3.4cm}{\color{mred}\raggedleft Plastic Caching\\\&
                                  Episodic-Like Memory\\\& Motivational Control}",
                 raw"\parbox{3cm}{\color{mblue}\raggedleft Episodic-Like Memory\\\&
                                  Motivational Control}",
                 raw"\color{mgreen} Motivational Control"],
                xticklabel_style = {rotate = 60, anchor = "east"}},
               [Plot(Coordinates(1:4, ys[:, i])) for i in 1:4]...)
pgfsave(joinpath(figpath, "logp_hat_rel_summary.tikz"), f2)


text:
# src/core.jl
# Core Definitions

import AbstractRationals: AbstractRational, missing_api, build_ratio
import Base: +, *, fma

export BoundType,
       isnonbinding,
       isbinding,
       mostbinding,
       AbstractBoundedRational,
       bound,
       isbounded,
       resolve

"""```
const BoundType
    === Union{<:Number, Missing, Nothing}
```
"""
const BoundType = Union{<:Number,Missing,Nothing}

"""```
isnonbinding(bound::BoundType)
```
"""
function isnonbinding(bound::BoundType)::Bool
    return isinf(bound) | ismissing(bound) | isnothing(bound) | isnan(bound)
end

"""```
isbinding(bound::BoundType)
```
"""
function isbinding(bound::BoundType)::Bool
    return !isinf(bound) & !ismissing(bound) & !isnothing(bound) & !isnan(bound)
end

"""```
mostbinding(bounds::BoundType...)
```
"""
function mostbinding(bounds::BoundType...)::BoundType
    components = filter(isbinding, bounds)
    return isempty(components) ? nothing : min(components)
end

"""```
abstract type AbstractBoundedRational{T, B <: BoundType}
    <: AbstractRational{T}
```
"""
abstract type AbstractBoundedRational{T,B<:BoundType} <: AbstractRational{T} end

"""```
bound(x::AbstractBoundedRational)
```
"""
function bound(x::AbstractBoundedRational{T,B})::B where {T,B}
    missing_api("bound", x)
end

"""```
isbounded(x::AbstractBoundedRational)
    === isbinding(bound(x))
```
"""
function isbounded(x::AbstractBoundedRational)
    return isbinding(bound(x))
end

"""```
resolve(num::T, den::T, bound::B)
```
"""
function resolve(num::T, den::T, bound::B) where {T<:Integer,B<:BoundType}
    if isnonbinding(bound)
        return num, den
    end
    scaled = num * bound
    n, r = divrem(scaled, den)
    return (r <= den - r ? n : n + one(n), bound)
end

"""```
x::AbstractBoundedRational * y::AbstractBoundedRational
```
"""
function *(x::AbstractBoundedRational, y::AbstractBoundedRational)
    return build_ratio(
        typeof(x),
        resolve(
            numerator(x) * numerator(y),
            denominator(x) * denominator(y),
            mostbinding(bound(x), bound(y)),
        )...,
    )
end

"""```
x::AbstractBoundedRational + y::AbstractBoundedRational
```
"""
function +(x::AbstractBoundedRational, y::AbstractBoundedRational)
    return build_ratio(
        typeof(x),
        resolve(
            numerator(x) * denominator(y) + numerator(y) * denominator(x),
            denominator(x) * denominator(y),
            mostbinding(bound(x), bound(y)),
        )...,
    )
end

"""```
fma(x::AbstractBoundedRational, y::AbstractBoundedRational, z::AbstractBoundedRational)
```
"""
function fma(x::AbstractBoundedRational, y::AbstractBoundedRational, z::AbstractBoundedRational)
    return build_ratio(
        typeof(x),
        resolve(
            numerator(x) * numerator(y),
            denominator(x) * denominator(y),
            mostbinding(bound(x), bound(y), bound(z)),
        )...,
    )
end


text:
# Stores data in a simple dict in memory
struct DictStore <: AbstractStore
  a::Dict{String,Vector{UInt8}}
end
DictStore() = DictStore(Dict{String,Vector{UInt8}}())

Base.show(io::IO,d::DictStore) = print(io,"Dictionary Storage")
_pdict(d::DictStore,p) = filter(((k,v),)->startswith(k,p),d.a)
function storagesize(d::DictStore,p)
  sum(i->last(split(i[1],'/')) ∉ (".zattrs",".zarray") ? sizeof(i[2]) : zero(sizeof(i[2])), _pdict(d,p))
end

function Base.getindex(d::DictStore,i::AbstractString)
  get(d.a,i,nothing)
end
function Base.setindex!(d::DictStore,v,i::AbstractString)
  d.a[i] = v
end
Base.delete!(d::DictStore, i::AbstractString) = delete!(d.a,i)

function subdirs(d::DictStore,p)
  d2 = _pdict(d,p)
  _searchsubdict(d2,p,(sp,lp)->length(sp) > lp+1)
end

function subkeys(d::DictStore,p)
  d2 = _pdict(d,p)
  _searchsubdict(d2,p,(sp,lp)->length(sp) == lp+1)
end

function _searchsubdict(d2,p,condition)
  o = Set{String}()
  pspl = split(rstrip(p,'/'),'/')
  lp = if length(pspl) == 1 && isempty(pspl[1])
    0
  else
    length(pspl)
  end
  for k in keys(d2)
    sp = split(k,'/')
    if condition(sp,lp)
      push!(o,sp[lp+1])
    end
  end
  collect(o)
end

#getsub(d::DictStore, p, n) = _substore(d,p).subdirs[n]

#path(d::DictStore) = ""


text:
using PGFPlots

using CSV
using DataFrames
using Missings
using Query

function only(a)
    @assert length(a) == 1
    return first(a)
end

# filename = Pkg.dir("Multilane", "data/baseline_curve_Saturday_24_Mar_21_08.csv")
# filename = Pkg.dir("Multilane", "data/baseline_curve_Sunday_25_Mar_12_09.csv")
# filename = Pkg.dir("Multilane", "data/baseline_curve_Wednesday_28_Mar_16_55.csv")
filename = Pkg.dir("Multilane", "data/baseline_curve_Friday_30_Mar_07_36.csv")

data = CSV.read(filename, nullable=true)

data[:time_to_lane] = data[:steps_to_lane].*data[:dt]
data[:brakes_per_step] = data[:nb_brakes]./data[:n_steps]

combined = by(data, [:solver, :lambda]) do df
    n = size(df, 1)
    return DataFrame(reached_lane=sum(df[:terminal].=="lane")/n,
                     had_brakes=sum(df[:nb_brakes].>=1)/n,
                    )
end

@show sort(combined, cols=[:solver])

plts = Plots.Plot[]
for df in groupby(combined, [:solver])
    df = sort(df, cols=[:lambda])
    name = only(unique(df[:solver]))
    plt = Plots.Linear(1.0.-df[:had_brakes], df[:reached_lane], legendentry=name)
    push!(plts, plt)
end
a = Axis(plts,
         ylabel="Fraction Successful",
         xlabel="Fraction Safe",
         legendPos="south west"
        )

file = tempname()*".pdf"
save(file, a)
run(`xdg-open $file`)


text:
# This file was generated, do not modify it. # hide
ensemble_model = EnsembleModel(model=knn_model, n=20);

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using OMEinsum, Test
using OMEinsum: SimpleBinaryRule, match_rule

@testset "analyse binary" begin
    size_dict = Dict(1=>1, 2=>2, 3=>3, 4=>4, 6=>6, 7=>7, 8=>8)
    c1, c2, cy, s1, s2, is, js, ys = OMEinsum.analyze_binary([1,2,3,4,8], [2,6,6,8,4,2], [7,2,1,2,2,6], size_dict)
    @test c1 == [1,4,8,2]
    @test c2 == [4,8,6,2]
    @test cy == [1,6,2]
    @test s1 == [1,32,2]
    @test s2 == [32,6,2]
    @test is == ['i', 'j', 'l']
    @test js == ['j', 'k', 'l']
    @test ys == ['i', 'k', 'l']
end

@testset "binary rules" begin
    size_dict = Dict(zip(('i', 'j', 'k', 'l'), ntuple(x->5, 4)))
    nmatch = 0
    for has_batch in [true, false]
        for i1 in [(), ('i',), ('j',), ('i','j'), ('j', 'i')]
            for i2 in [(), ('k',), ('j',), ('k','j'), ('j', 'k')]
                for i3 in [(), ('i',), ('k',), ('i','k'), ('k', 'i')]
                    @info i1, i2, i3, has_batch
                    i1_ = has_batch ? (i1..., 'l') : i1
                    i2_ = has_batch ? (i2..., 'l') : i2
                    i3_ = has_batch ? (i3..., 'l') : i3
                    a = randn(fill(5, length(i1_))...) |> OMEinsum.asarray
                    b = randn(fill(5, length(i2_))...) |> OMEinsum.asarray
                    code = EinCode((i1_,i2_),i3_)
                    rule = match_rule(code)
                    if rule isa SimpleBinaryRule
                        nmatch += 1
                        @test einsum(rule, (a, b)) ≈ loop_einsum(code, (a, b), size_dict)
                    else
                        @test einsum(code, (a, b), size_dict) ≈ loop_einsum(code, (a, b), size_dict)
                    end
                end
            end
        end
    end
    @test nmatch == 36
end

@testset "match binary rules" begin
    for code in [ein",->", ein",k->k", ein"i,->i", ein"j,j->", ein"i,k->ik", ein"i,k->ki",
                ein"j,jk->k", ein"j,kj->k", ein"ji,j->i", ein"ij,j->i", ein"ji,jk->ik",
                ein"ji,kj->ik", ein"ji,jk->ki", ein"ji,kj->ki", ein"ij,jk->ik", ein"ij,kj->ik", ein"ij,jk->ki", ein"ij,kj->ki"]
        @test OMEinsum.match_rule(code) == SimpleBinaryRule(code)
    end
    OMEinsum.match_rule(ein"ab,bc->ac") == SimpleBinaryRule(('i', 'j'), ('j', 'k'), ('i', 'k'))
    OMEinsum.match_rule(ein"ab,bce->ac") == OMEinsum.DefaultRule()
end

@testset "binary einsum" begin
    for code in [
        ein",->",
        ein"ijl,jl->il",
        ein"ab,bc->ac",
        ein"abb,bc->ac",  # with diag in
        ein"ab,bc->acc",  # with diag out
        ein"ab,bce->ac",  # with sum in
        ein"ab,bc->ace",  # with sum out
        ein"bal,bcl->lcae",  # with perm in
        ein"bal,bcl->ca",  # with multi-edge in
        ein"bal,bc->lca",  # with multi-edge out
        ein"ddebal,bcf->lcac",  # with all
    ]
        xs = [OMEinsum.asarray(randn(ComplexF64, fill(5, length(ix))...)) for ix in OMEinsum.getixs(code)]
        size_dict = Dict(zip(('a', 'b', 'c', 'd', 'e', 'f','i','j','k','l'), ntuple(x->5, 10)))
        @test einsum(code, (xs...,), size_dict) ≈ loop_einsum(code, (xs...,), size_dict)
    end
end

@testset "regression test" begin
    x, y = randn(fill(2, 8)...), randn(fill(2, 5)...)
    code = EinCode(((10, 5, 4, 2, 6, 6, 7, 9), (1, 9, 8, 2, 10)), (4, 10, 9, 7, 10))
    @test code(x, y) ≈ loop_einsum(code, (x, y), OMEinsum.get_size_dict(OMEinsum.getixs(code), (x, y)))
end


text:
module CompileUtils

using ..AExpressions
using Distributions: Categorical
using MLStyle: @match

export AutumnCompileError, compile, compilestatestruct, compileinitstate, compileinitnext, compileprevfuncs, compilebuiltin

"Autumn Compile Error"
struct AutumnCompileError <: Exception
  msg
end
AutumnCompileError() = AutumnCompileError("")

# ----- Compile Helper Functions ----- #

function compile(expr::AExpr, data::Dict{String, Any}, parent::Union{AExpr, Nothing}=nothing)
  arr = [expr.head, expr.args...]
  res = @match arr begin
    [:if, args...] => :($(compile(args[1], data)) ? $(compile(args[2], data)) : $(compile(args[3], data)))
    [:assign, args...] => compileassign(expr, data, parent)
    [:typedecl, args...] => compiletypedecl(expr, data, parent)
    [:external, args...] => compileexternal(expr, data)
    [:let, args...] => compilelet(expr, data)
    [:case, args...] => compilecase(expr, data)
    [:typealias, args...] => compiletypealias(expr, data)
    [:lambda, args...] => :($(compile(args[1], data)) -> $(compile(args[2], data)))
    [:list, args...] => :([$(map(x -> compile(x, data), expr.args)...)])
    [:call, args...] => compilecall(expr, data)
    [:field, args...] => :($(compile(expr.args[1], data)).$(compile(expr.args[2], data)))
    [args...] => throw(AutumnCompileError(string("Invalid AExpr Head: ", expr.head))) # if expr head is not one of the above, throw error
  end
end

function compile(expr::AbstractArray, data::Dict{String, Any}, parent::Union{AExpr, Nothing}=nothing)
  if length(expr) == 0 || (length(expr) > 1 && expr[1] != :List)
    throw(AutumnCompileError("Invalid List Syntax"))
  elseif expr[1] == :List
    :(Array{$(compile(expr[2:end], data))})
  else
    expr[1]
  end
end

function compile(expr, data::Dict{String, Any}, parent::Union{AExpr, Nothing}=nothing)
  expr
end

function compileassign(expr::AExpr, data::Dict{String, Any}, parent::Union{AExpr, Nothing})
  # get type, if declared
  type = haskey(data["types"], expr.args[1]) ? data["types"][expr.args[1]] : nothing
  if (typeof(expr.args[2]) == AExpr && expr.args[2].head == :fn)
    if type !== nothing # handle function with typed arguments/return type
      args = compile(expr.args[2].args[1], data).args # function args
      argtypes = map(x -> compile(x, data), type.args[1:(end-1)]) # function arg types
      tuples = [(args[i], argtypes[i]) for i in [1:length(args);]]
      typedargexprs = map(x -> :($(x[1])::$(x[2])), tuples)
      quote
        function $(compile(expr.args[1], data))($(typedargexprs...))::$(compile(type.args[end], data))
          $(compile(expr.args[2].args[2], data))
        end
      end
    else # handle function without typed arguments/return type
      quote
        function $(compile(expr.args[1], data))($(compile(expr.args[2].args[1], data).args[2]...))
            $(compile(expr.args[2].args[2], data))
        end
      end
    end
  else # handle non-function assignments
    # handle global assignments
    if parent !== nothing && (parent.head == :program)
      if (typeof(expr.args[2]) == AExpr && expr.args[2].head == :initnext)
        push!(data["initnext"], expr)
      else
        push!(data["lifted"], expr)
      end
      :()
    # handle non-global assignments
    else
      if type !== nothing
        :($(compile(expr.args[1], data))::$(compile(type, data)) = compile(expr.args[2], data))
      else
          :($(compile(expr.args[1], data)) = $(compile(expr.args[2], data)))
      end
    end
  end
end

function compiletypedecl(expr::AExpr, data::Dict{String, Any}, parent::Union{AExpr, Nothing})
  if (parent !== nothing && (parent.head == :program || parent.head == :external))
    data["types"][expr.args[1]] = expr.args[2]
    :()
  else
    :(local $(compile(expr.args[1], data))::$(compile(expr.args[2], data)))
  end
end

function compileexternal(expr::AExpr, data::Dict{String, Any})
  push!(data["external"], expr.args[1])
  compiletypedecl(expr.args[1], data, expr)
end

function compiletypealias(expr::AExpr, data::Dict{String, Any})
  name = expr.args[1]
  fields = map(field -> (
    :($(field.args[1])::$(field.args[2]))
  ), expr.args[2].args)
  quote
    struct $(name)
      $(fields...)
    end
  end
end

function compilecall(expr::AExpr, data::Dict{String, Any})
  fnName = expr.args[1]
  if !(fnName in binaryOperators) && fnName != :prev
    :($(fnName)($(map(x -> compile(x, data), expr.args[2:end])...)))
  elseif fnName == :prev
    :($(Symbol(string(expr.args[2]) * "Prev"))($(map(compile, expr.args[3:end])...)))
  elseif fnName != :(==)
    :($(fnName)($(compile(expr.args[2], data)), $(compile(expr.args[3], data))))
  else
    :($(compile(expr.args[2], data)) == $(compile(expr.args[3], data)))
  end
end

function compilelet(expr::AExpr, data::Dict{String, Any})
  quote
    $(map(x -> compile(x, data), expr.args)...)
  end
end

function compilecase(expr::AExpr, data::Dict{String, Any})
  quote
    @match $(compile(expr.args[1], data)) begin
      $(map(x -> :($(compile(x.args[1], data)) => $(compile(x.args[2], data))), expr.args[2:end])...)
    end
  end
end

function compileinitnext(data::Dict{String, Any})
  initFunction = quote
    function init($(map(x -> :($(compile(x.args[1], data))::Union{$(compile(data["types"][x.args[1]], data)), Nothing}), data["external"])...))::STATE
      $(compileinitstate(data))
      $(map(x -> :($(compile(x.args[1], data)) = $(compile(x.args[2].args[1], data))), data["initnext"])...)
      $(map(x -> :($(compile(x.args[1], data)) = $(compile(x.args[2], data))), data["lifted"])...)
      $(map(x -> :(state.$(Symbol(string(x.args[1])*"History"))[state.time] = $(compile(x.args[1], data))),
            vcat(data["external"], data["initnext"], data["lifted"]))...)
      deepcopy(state)
    end
    end
  nextFunction = quote
    function next($([:(old_state::STATE), map(x -> :($(compile(x.args[1], data))::Union{$(compile(data["types"][x.args[1]], data)), Nothing}), data["external"])...]...))::STATE
      global state = deepcopy(old_state)
      state.time = state.time + 1
      $(map(x -> :($(compile(x.args[1], data)) = $(compile(x.args[2].args[2], data))), data["initnext"])...)
      $(map(x -> :($(compile(x.args[1], data)) = $(compile(x.args[2], data))), data["lifted"])...)
      $(map(x -> :(state.$(Symbol(string(x.args[1])*"History"))[state.time] = $(compile(x.args[1], data))),
            vcat(data["external"], data["initnext"], data["lifted"]))...)
      deepcopy(state)
    end
    end
    [initFunction, nextFunction]
end

# construct STATE struct
function compilestatestruct(data::Dict{String, Any})
  stateParamsInternal = map(expr -> :($(Symbol(string(expr.args[1]) * "History"))::Dict{Int64, $(haskey(data["types"], expr.args[1]) ? compile(data["types"][expr.args[1]], data) : Any)}),
                            vcat(data["initnext"], data["lifted"]))
  stateParamsExternal = map(expr -> :($(Symbol(string(expr.args[1]) * "History"))::Dict{Int64, Union{$(compile(data["types"][expr.args[1]], data)), Nothing}}),
                            data["external"])
  quote
    mutable struct STATE
      time::Int
      $(stateParamsInternal...)
      $(stateParamsExternal...)
    end
  end
end

# initialize state::STATE variable
function compileinitstate(data::Dict{String, Any})
  initStateParamsInternal = map(expr -> :(Dict{Int64, $(haskey(data["types"], expr.args[1]) ? compile(data["types"][expr.args[1]], data) : Any)}()),
                                vcat(data["initnext"], data["lifted"]))
  initStateParamsExternal = map(expr -> :(Dict{Int64, Union{$(compile(data["types"][expr.args[1]], data)), Nothing}}()),
                                data["external"])
  initStateParams = [0, initStateParamsInternal..., initStateParamsExternal...]
  initState = :(state = STATE($(initStateParams...)))
  initState
end

# ----- Built-In and Prev Function Helpers ----- #

function compileprevfuncs(data::Dict{String, Any})
  prevFunctions = map(x -> quote
        function $(Symbol(string(x.args[1]) * "Prev"))(n::Int=1)::$(haskey(data["types"], x.args[1]) ? compile(data["types"][x.args[1]], data) : Any)
          state.$(Symbol(string(x.args[1]) * "History"))[state.time - n >= 0 ? state.time - n : 0]
        end
        end,
  vcat(data["initnext"], data["lifted"]))
  prevFunctions = vcat(prevFunctions, map(x -> quote
        function $(Symbol(string(x.args[1]) * "Prev"))(n::Int=1)::Union{$(compile(data["types"][x.args[1]], data)), Nothing}
          state.$(Symbol(string(x.args[1]) * "History"))[state.time - n >= 0 ? state.time - n : 0]
        end
        end,
  data["external"]))
  prevFunctions
end

function compilebuiltin()
  occurredFunction = builtInDict["occurred"]
  uniformChoiceFunction = builtInDict["uniformChoice"]
  uniformChoiceFunction2 = builtInDict["uniformChoice2"]
  minFunction = builtInDict["min"]
  clickType = builtInDict["clickType"]
  rangeFunction = builtInDict["range"]
  [occurredFunction, uniformChoiceFunction, uniformChoiceFunction2, minFunction, clickType, rangeFunction]
end

const builtInDict = Dict([
"occurred"        =>  quote
                        function occurred(click)
                          click !== nothing
                        end
                      end,
"uniformChoice"   =>  quote
                        function uniformChoice(freePositions)
                          freePositions[rand(Categorical(ones(length(freePositions))/length(freePositions)))]
                        end
                      end,
"uniformChoice2"   =>  quote
                        function uniformChoice(freePositions, n)
                          map(idx -> freePositions[idx], rand(Categorical(ones(length(freePositions))/length(freePositions)), n))
                        end
                      end,
"min"              => quote
                        function min(arr)
                          min(arr...)
                        end
                      end,
"clickType"       =>  quote
                        struct Click
                          x::BigInt
                          y::BigInt
                        end
                      end,
"range"           => quote
                      function range(start::BigInt, stop::BigInt)
                        [start:stop;]
                      end
                     end
])

# binary operators
const binaryOperators = [:+, :-, :/, :*, :&, :|, :>=, :<=, :>, :<, :(==), :!=, :%, :&&]

end


text:
"""
    envelope_lattice(model::MetabolicModel, rids::Vector{String}; kwargs...)

Version of [`envelope_lattice`](@ref) that works on string reaction IDs instead
of integer indexes.
"""
envelope_lattice(model::MetabolicModel, rids::Vector{String}; kwargs...) =
    envelope_lattice(model, Vector{Int}(indexin(rids, reactions(model))); kwargs...)

"""
    envelope_lattice(
        model::MetabolicModel,
        ridxs::Vector{Int};
        samples = 10,
        ranges = collect(zip(bounds(model)...))[ridxs],
        reaction_samples = fill(samples, length(ridxs)),
    )

Create a lattice (list of "tick" vectors) for reactions at indexes `ridxs` in a
model. Arguments `samples`, `ranges`, and `reaction_samples` may be optionally
specified to customize the lattice creation process.
"""
envelope_lattice(
    model::MetabolicModel,
    ridxs::Vector{Int};
    samples = 10,
    ranges = collect(zip(bounds(model)...))[ridxs],
    reaction_samples = fill(samples, length(ridxs)),
) = (
    lb .+ (ub - lb) .* ((1:s) .- 1) ./ max(s - 1, 1) for
    (s, (lb, ub)) in zip(reaction_samples, ranges)
)

"""
    objective_envelope(model::MetabolicModel, rids::Vector{String}, args...; kwargs...)

Version of [`objective_envelope`](@ref) that works on string reaction IDs
instead of integer indexes.
"""
objective_envelope(model::MetabolicModel, rids::Vector{String}, args...; kwargs...) =
    objective_envelope(
        model,
        Vector{Int}(indexin(rids, reactions(model))),
        args...;
        kwargs...,
    )

"""
    objective_envelope(
        model::MetabolicModel,
        ridxs::Vector{Int},
        optimizer;
        modifications = [],
        lattice_args = (),
        lattice = envelope_lattice(model, ridxs; lattice_args...),
        analysis = screen_optimize_objective,
        kwargs...,
    )

Compute an array of objective values for the `model` for rates of reactions
specified `ridxs` fixed to a regular range of values between their respective
lower and upper bounds.

This can be used to compute a "production envelope" of a metabolite; but
generalizes to any specifiable objective and to multiple dimensions of the
examined space. For example, to retrieve a production envelope of any
metabolite, set the objective coefficient vector of the `model` to a vector
that contains a single `1` for the exchange reaction that "outputs" this
metabolite, and run [`objective_envelope`](@ref) with the exchange reaction of
the "parameter" metabolite specified in `ridxs`.

Returns a named tuple that contains `lattice` with reference values of the
metabolites, and an N-dimensional array `values` with the computed objective
values, where N is the number of specified reactions.  Because of the
increasing dimensionality, the computation gets very voluminous with increasing
length of `ridxs`. The `lattice` for computing the optima can be supplied in
the argument; by default it is created by [`envelope_lattice`](@ref) called on
the model and reaction indexes. Additional arguments for the call to
[`envelope_lattice`](@ref) can be optionally specified in `lattice_args`.

`kwargs` are internally forwarded to [`screen_optmodel_modifications`](@ref).
`modifications` are appended to the list of modifications after the lattice
bounds are set. By default, this returns the objective values for all points in
the lattice; alternate outputs can be implemented via the `analysis` argument.

# Example
```
julia> m = load_model("test/downloaded/e_coli_core.xml");

julia> envelope = objective_envelope(m, ["R_EX_gln__L_e", "R_EX_fum_e"],
                                     Tulip.Optimizer;
                                     lattice_args=(samples=6,));

julia> envelope.lattice   # the reaction rates for which the optima were computed
2-element Vector{Vector{Float64}}:
 [0.0, 200.0, 400.0, 600.0, 800.0, 1000.0]
 [0.0, 200.0, 400.0, 600.0, 800.0, 1000.0]

julia> envelope.values   # the computed flux objective values for each reaction rate combination
6×6 Matrix{Float64}:
  0.873922   9.25815  17.4538  19.56   20.4121  20.4121
 13.0354    17.508    19.9369  21.894  22.6825  22.6825
 16.6666    18.6097   20.2847  21.894  22.6825  22.6825
 16.6666    18.6097   20.2847  21.894  22.6825  22.6825
 16.6666    18.6097   20.2847  21.894  22.6825  22.6825
 16.6666    18.6097   20.2847  21.894  22.6825  22.6825
```
"""
objective_envelope(
    model::MetabolicModel,
    ridxs::Vector{Int},
    optimizer;
    modifications = [],
    lattice_args = (),
    lattice = envelope_lattice(model, ridxs; lattice_args...),
    analysis = screen_optimize_objective,
    kwargs...,
) = (
    lattice = collect.(lattice),
    values = screen_optmodel_modifications(
        model,
        optimizer;
        modifications = collect(
            vcat(
                [
                    (_, optmodel) -> begin
                        for (i, ridx) in enumerate(ridxs)
                            set_normalized_rhs(optmodel[:lbs][ridx], fluxes[i])
                            set_normalized_rhs(optmodel[:ubs][ridx], fluxes[i])
                        end
                    end,
                ],
                modifications,
            ) for fluxes in Iterators.product(lattice...)
        ),
        analysis = analysis,
        kwargs...,
    ),
)


text:
include("../src/solution.jl")
include("../src/solvers/explore.jl")

# ===========
# Préparartion des arguments
args = Args.parse_commandline("explore")
const TEST_LEVEL = 1
Args.set(:level, TEST_LEVEL) # sera passé temporairement à 0
# Args.set(:itermax, 1_000) # Inutile pour ExploreSolver
# itermax = Args.get(:itermax)

# ===========
# Construction de l'instance mini
inst = instance_build_mini10()
@test inst.nb_planes == 10

# ===========
lg1("1. Création d'un ExploreSolver... ")
# Args.set(:level, 0)
# Args.show_args()
solver = ExploreSolver(inst)
# println(Args.to_s_dict(Args.args))
# shuffle!(solver.cursol)
ln1("fait.")

# ===========
lg1("2. Mélange de la solution initiale... ")
initial_sort!(solver.cursol, presort=:shuffle)
copy!(solver.bestsol, solver.cursol)
ln1(" => (initcost=$(solver.cursol.cost)) fait")

# ===========
itermax = 10_000
lg1("3. Résolution avec itermax=$itermax) itérations ")
lg1("(cost=700 <= 5500.0 ?)... ")
# solve(solver, Args.get(:itermax))
Args.set(:level, 0)
solve(solver, itermax, small_shift!)
Args.set(:level, TEST_LEVEL)
bestsol = solver.bestsol
@test bestsol.cost <= 5500.0
ln1(" => ok ($(bestsol.cost))")

# ===========
itermax = 100_000
lg1("4. Résolution avec itermax=$itermax) itérations supplémentaires")
lg1("(cost=700 <= 5500.0 ?)... ")
# solve(solver, Args.get(:itermax))
Args.set(:level, 0)
solve(solver, itermax, swap!)
Args.set(:level, TEST_LEVEL)
bestsol = solver.bestsol
@test bestsol.cost <= 5500.0
ln1(" => ok ($(bestsol.cost))\n")

#./