text:
#=
    The consensus tree functions are implemented on basis of the following paper:

   Jesper Jansson, Chuanqi Shen, and Wing-Kin Sung. 2016. Improved algorithms
   for constructing consensustrees. J. ACM 63, 3, Article 28 (June 2016), 24 pages

   Link: https://dl.acm.org/doi/pdf/10.1145/2925985

   Most of the major functions include a note to find the section of the paper
   they were based on.
=#

"""
    majority_consensus_tree(trees::Vector{T}, percentage::Float64=0.5)
        ::T where T<:AbstractNode

Construct the majority rule consensus tree from a set of trees that share the
same leafset. By default the output tree includes clusters that occur in over
50% of the trees. This can be customized when calling the function. The function
returns the root node of the majority consensus tree, from which it can be
traversed. The algorithm is based on section 3 and 6.1 of:

Jesper Jansson, Chuanqi Shen, and Wing-Kin Sung. 2016. Improved algorithms
for constructing consensustrees. J. ACM 63, 3, Article 28 (June 2016), 24 pages
https://dl.acm.org/doi/pdf/10.1145/2925985
"""
function majority_consensus_tree(trees::Vector{T}, percentage::Float64=0.5)::T where T<:AbstractNode
    merged_tree = deepcopy(trees[1])
    check_leafsets(trees)
    nodes = post_order(merged_tree)
    # save leaf ranks to order the resulting tree in the end
    leaf_ranks = get_leaf_ranks(nodes)
    node_counts = convert(Vector{Int64}, ones(length(nodes)))
    count_dict = Dict(zip(nodes, node_counts))
    for tree in trees[2:end]
        nodes = level_order(merged_tree)
        is_common_cluster = find_common_clusters(tree, merged_tree)
        for node in nodes
            # increment count of clusters of the merged_tree that are in the other tree
            if is_common_cluster[node.num][1] == true
                count_dict[node] += 1
            # delete clusters which are not
            else
                count_dict[node] -= 1
                if count_dict[node] == 0
                    delete_node!(node)
                end # if
            end # else
        end # for
        set_binary!(merged_tree)
        number_nodes!(merged_tree)
        compatible_tree = one_way_compatible(tree, merged_tree)
        inserted_nodes = merge_trees!(compatible_tree, merged_tree)
        # intialize counts for the new nodes
        for node in inserted_nodes
            count_dict[node] = 1
        end # for
    end # for
    # find clusters of the final tree in the other tree, store inc_length and count
    set_node_stats!(merged_tree, trees, true, percentage)
    # order the resulting tree
    nodes = post_order(merged_tree)
    cluster_start_indeces = Dict{Node, Int64}()
    for node in nodes
        if node.nchild != 0
            cluster_start_indeces[node] = cluster_start_indeces[node.children[1]]
        else
            cluster_start_indeces[node] = leaf_ranks[node.name]
        end # if/else
    end # for
    order_tree!(merged_tree, cluster_start_indeces)
    return merged_tree
end

"""
    loose_consensus_tree(trees::Vector{T})::T where T<:AbstractNode

Construct the loose consensus tree from a set of trees that share the same
leafset. I.e. a tree with all the clusters that appear in at least one tree
and are compatible with all trees. Returns the root node of the loose consensus
tree, from which it can be traversed. This algorithm is based on section 4 and
6.1 of:

Jesper Jansson, Chuanqi Shen, and Wing-Kin Sung. 2016. Improved algorithms
for constructing consensustrees. J. ACM 63, 3, Article 28 (June 2016), 24 pages
https://dl.acm.org/doi/pdf/10.1145/2925985
"""
function loose_consensus_tree(trees::Vector{T})::T where T<:AbstractNode
    r_tree = trees[1]
    check_leafsets(trees)
    nodes = post_order(r_tree)
    # save leaf ranks to order the resulting tree in the end
    leaf_ranks = get_leaf_ranks(nodes)
    for i in 2:length(trees)
        compatible_tree = one_way_compatible(r_tree, trees[i])
        r_tree = deepcopy(trees[i])
        merge_trees!(compatible_tree, r_tree)
    end
    for tree in trees
        r_tree = one_way_compatible(r_tree, tree)
    end
    set_node_stats!(r_tree, trees, false)
    # order the resulting tree
    nodes = post_order(r_tree)
    cluster_start_indeces = Dict{Node, Int64}()
    for node in nodes
        if node.nchild != 0
            cluster_start_indeces[node] = cluster_start_indeces[node.children[1]]
        else
            cluster_start_indeces[node] = leaf_ranks[node.name]
        end # if/else
    end # for
    order_tree!(r_tree, cluster_start_indeces)
    return r_tree
end

"""
    greedy_consensus_tree(trees::Vector{T})::T where T<:AbstractNode

Construct the greedy consensus tree from a set of trees that share the same
leafset.  Returns the root node of the greedy consensus tree, from which it can
be traversed. This algorithm is based on section 5 and 6.1 of:

Jesper Jansson, Chuanqi Shen, and Wing-Kin Sung. 2016. Improved algorithms
for constructing consensustrees. J. ACM 63, 3, Article 28 (June 2016), 24 pages
https://dl.acm.org/doi/pdf/10.1145/2925985
"""
function greedy_consensus_tree(trees::Vector{T})::T where T<:AbstractNode
    leafset = Set([n.name for n in get_leaves(trees[1])])
    l = length(leafset)
    check_leafsets(trees)
    nodes = post_order(trees[1])
    # PSEUDOCODE Step 1
    leaf_ranks = get_leaf_ranks(nodes)
    bit_vectors = Vector{BitVector}()
    for tree in trees
        temp_bit_vectors = Dict{Node, BitVector}()
        for node in post_order(tree)
            if node.root
                continue
            elseif node.nchild == 0
                bit_vector = falses(l)
                bit_vector[leaf_ranks[node.name]] = 1
                temp_bit_vectors[node] = bit_vector
            else
                bit_vector = temp_bit_vectors[node.children[1]]
                for child in node.children[2:end]
                    bit_vector = bit_vector .| temp_bit_vectors[child]
                end # for
                temp_bit_vectors[node] = bit_vector
            end # if/else
        end # for
        append!(bit_vectors, values(temp_bit_vectors))
    end # for
    # PSEUDOCODE Step 2
    sort!(bit_vectors, alg=QuickSort)
    # PSEUDOCODE Step 3
    cluster_queue = count_cluster_occurences(bit_vectors)
    # PSEUDOCODE Step 4
    greedy_consensus_tree = Node()
    for leaf in get_leaves(trees[1])
        add_child!(greedy_consensus_tree, deepcopy(leaf))
    end # for
    # PSEUDOCODE Step 5
    while !isempty(cluster_queue)
        cluster = dequeue!(cluster_queue)
        # ignore trivial leaf clusters
        sum(cluster) == 1 && continue
        cluster_length = sum(cluster)
        # num stores the number of shared leaves between a node and the current
        # cluster. It also stores the number of leaf nodes of the node.
        num = Dict{Node, Vector{Int64}}()
        for node in post_order(greedy_consensus_tree)
            if node.nchild == 0
                cluster[leaf_ranks[node.name]] == 1 ? num[node] = [1, 1] :
                                                      num[node] = [0, 1]
            else
                num[node] = [sum([num[child][1] for child in node.children])]
                n_leaves = 0
                for child in node.children
                        n_leaves += num[child][2]
                end # for
                push!(num[node], n_leaves)
                children = Vector{Node}()
                if num[node][1] == cluster_length
                    insert = true
                    for child in node.children
                        if num[child][1] == num[child][2]
                            push!(children, child)
                        elseif num[child][1] != 0
                            insert = false
                            break
                        end # if
                    end # for
                    insert && insert_node!(node, children)
                    number_nodes!(greedy_consensus_tree)
                    break
                end # if
            end # if/else
        end # for
    end # while
    set_node_stats!(greedy_consensus_tree, trees, false)
    return greedy_consensus_tree
end # greedy_consensus_tree

"""
    count_cluster_occurences(bit_vectors::BitVector)
        ::PriorityQueue{BitVector, Int64}

--- INTERNAL ---
Helper function for greedy_consensus_tree that counts the occurrences of each
bit vector representing a cluster. Returns a priority queue.
"""
function count_cluster_occurences(bit_vectors::Vector{BitVector}
                                 )::PriorityQueue{BitVector, Int64}

    cluster_queue = PriorityQueue{BitVector, Int64}()
    start = bit_vectors[1]
    count = 1
    for bit_vector in bit_vectors[2:end]
        if bit_vector == start
            count += 1
        else
            cluster_queue[start] = -count
            start = bit_vector
            count = 1
        end # if/else
    end # for
    cluster_queue[start] = -count
    return cluster_queue
end # count_cluster_occurences

"""
    set_node_stats!(tree::T, trees::Vector{T}, majority::Bool)
        ::Nothing where T<:AbstractNode

--- INTERNAL ---
Helper function for the construction of a consensus tree. Calculates the
inc_lengths and statistics of the nodes in the consensus tree. If dealing with a
(in-progress) majority consensus tree, this function will also delete its
non-majority clusters (handled by the 'majority' boolean).
"""
function set_node_stats!(main_tree::T, trees::Vector{T}, majority::Bool,
                         percentage::Float64=0.0)::Nothing where T<:AbstractNode

    nodes = level_order(main_tree)
    count_dict = Dict(zip([n.num for n in nodes], zeros(Int64, length(nodes))))
    inc_length_dict = Dict{Int64, Vector{Float64}}()
    for tree in trees
        is_common_cluster = find_common_clusters(tree, main_tree)
        for node in nodes
            if is_common_cluster[node.num][1]
                count_dict[node.num] += 1
                try
                    push!(inc_length_dict[node.num], is_common_cluster[node.num][2])
                catch KeyError
                    inc_length_dict[node.num] = [is_common_cluster[node.num][2]]
                end # try/catch
            end # if
        end # for
    end # for
    half = length(trees) * percentage
    # delete non-majority clusters
    for node in nodes
        if majority && count_dict[node.num] <= half
            delete_node!(node)
        else
            node.inc_length = mean(inc_length_dict[node.num])
            node.stats["sd"] = std(inc_length_dict[node.num])
            node.stats["median"] = median(inc_length_dict[node.num])
            node.stats["frequency"] = count_dict[node.num] / length(trees)
        end # if
    end # for
end # set_node_stats

"""
    find_common_clusters(ref_tree, tree:T)
        ::Dict{Int64, Tuple{Bool, Union{Float64, Missing}}}

--- INTERNAL ---
Use Day's algorithm to create a dictionary, that tells us for each node of the
second input tree, if its corresponding cluster is a common cluster of the trees.

Returns dictionary.

* `ref_tree` : Tree used to create dictionary.

* `tree` : Tree used to create dictionary.

"""
function find_common_clusters(ref_tree::T, tree::T) where T<:AbstractNode
    ref_nodes = post_order(ref_tree)
    leaves_dict = Dict{String, Tuple{Int64, Float64}}()
    clusters = Dict{Node, Vector{String}}()
    cluster_dict = Dict{Tuple{Int64, Int64}, Int64}()
    dict_count = 1
    for ref_node in ref_nodes
        if ref_node.nchild == 0
            leaves_dict[ref_node.name] = (dict_count, ref_node.inc_length)
            dict_count += 1
        else
            cluster = Vector{String}()
            for child in ref_node.children
                child.nchild == 0 ? push!(cluster, child.name) :
                                    append!(cluster, clusters[child])
            end # for
            clusters[ref_node] = cluster
            if ref_node.root == true || ref_node == ref_node.mother.children[1]
                last_index = leaves_dict[last(cluster)][1]
                first_index = leaves_dict[first(cluster)][1]
                cluster_dict[(first_index, last_index)] = last_index
            else
                last_index = leaves_dict[last(cluster)][1]
                first_index = leaves_dict[first(cluster)][1]
                cluster_dict[(first_index, last_index)] = first_index
            end # if/else
        end # if/else
    end # for
    for value in values(clusters)
        sort!(value)
    end # for
    clusters_reverse = Dict(value => key.inc_length for (key, value) in clusters)
    is_common_cluster = Dict{Int64, Tuple{Bool, Union{Float64, Missing}}}()
    clusters = Dict{Node, Vector{String}}()
    leaf_count = 0
    nodes = post_order(tree)
    for node in nodes
        if node.nchild == 0
            leaf_count += 1
            try
                leaves_dict[node.name]
            catch KeyError
                throw(ArgumentError(
                      "The leafs sets of the trees need to be identical"
                      ))
            end # try
            is_common_cluster[node.num] = (true, leaves_dict[node.name][2])
        else
            cluster = Vector{String}()
            for child in node.children
                child.nchild == 0 ? push!(cluster, child.name) :
                                    append!(cluster, clusters[child])
            end # for
            clusters[node] = cluster
            cluster_indeces = [leaves_dict[leaf][1] for leaf in cluster]
            if length(cluster) != maximum(cluster_indeces) - minimum(cluster_indeces) + 1
                is_common_cluster[node.num] = (false, missing)
            elseif haskey(cluster_dict, (minimum(cluster_indeces), maximum(cluster_indeces)))
                is_common_cluster[node.num] = (true, clusters_reverse[sort!(cluster)])
            else
                is_common_cluster[node.num] = (false, missing)
            end # if/else
        end # if/else
    end # for
    if leaf_count != length(keys(leaves_dict))
        throw(ArgumentError("The leafs sets of the trees need to be identical"))
    end # if
    return is_common_cluster
end # find_common_clusters

"""
    one_way_compatible(ref_tree::T, tree::T)::T where T<:AbstractNode

--- INTERNAL ---
Takes two trees and returns a copy of the first one, where all the clusters that
are not compatible with the second tree are removed.

Returns root node of edited tree.

* `ref_tree` : tree used to determine clusters.

* `tree` : tree used to determine clusters.
"""
function one_way_compatible(ref_tree::T, tree::T)::T where T<:AbstractNode
    ref_tree_copy = deepcopy(ref_tree)
    ref_nodes = post_order(ref_tree_copy)
    cluster_start_indeces = get_cluster_start_indeces(ref_nodes, tree)
    leaves::Vector{Node} = order_tree!(tree, cluster_start_indeces)
    leaf_ranks_reverse = Dict(node.name => ind for (ind, node) in enumerate(leaves))
    xleft_dict, xright_dict = depth_dicts(leaves)
    marked_nodes = Dict{Int64, Bool}()
    for ref_node in ref_nodes
        check_node!(ref_node, leaves, leaf_ranks_reverse, xleft_dict,
                    xright_dict, marked_nodes)
    end # for
    for node in level_order(ref_tree_copy)
        if node.nchild != 0 && marked_nodes[node.num] == true
            delete_node!(node)
        end # if
    end # for
    set_binary!(ref_tree_copy)
    number_nodes!(ref_tree_copy)
    return ref_tree_copy
end # function one_way_compatible

"""
    merge_trees!(ref_tree::T, tree::T)::Vector{T}} where T<:AbstractNode

--- INTERNAL ---
Merge two compatible trees, i.e. inserts all cluster of the first tree, which
aren't already in the second tree, into the second tree

Returns a Tuple containing the root node of the merged tree, as well as a vector
of nodes that were added.

* `ref_tree` : tree used to determine common clusters.

* `tree` : tree used to determine common clusters.
"""
function merge_trees!(ref_tree::T, tree::T)::Vector{T} where T<:AbstractNode
    ref_nodes = post_order(ref_tree)
    cluster_start_indeces = get_cluster_start_indeces(ref_nodes, tree)
    leaves::Vector{Node} = order_tree!(tree, cluster_start_indeces)
    leaf_ranks_reverse = Dict(node.name => ind for (ind, node) in enumerate(leaves))
    xleft_dict, xright_dict = depth_dicts(leaves)
    inserted_nodes = Vector{Node}()
    count = -1
    for ref_node in ref_nodes
        tuple = check_node!(ref_node, leaves, leaf_ranks_reverse, xleft_dict,
                            xright_dict)
        isnothing(tuple) && continue
        d, e, r, start_node, stop_node, left_path, right_path = tuple
        left = d == r.children[1]
        right = e == r.children[end]
        if !(left && right)
            index_d = findfirst(x -> x == d, r.children)
            index_e = findfirst(x -> x == e, r.children)
            inserted_node = insert_node!(r, r.children[index_d:index_e])
            push!(inserted_nodes, inserted_node)
            # ensures correct depth
            inserted_node.binary = string(r.binary, ",z")
            # give unique number to avoid false positive "==" statements
            inserted_node.num = count
            count -= 1
            inserted_depth = node_depth(inserted_node)
            if !left
                left_path[inserted_depth] = inserted_node
                right_path[inserted_depth] = inserted_node
                xleft_dict[start_node] = (inserted_node, left_path)
            end
            if !right
                right_path[inserted_depth] = inserted_node
                left_path[inserted_depth] = inserted_node
                xright_dict[stop_node] = (inserted_node, right_path)
            end # if
        end # if
    end # for
    set_binary!(tree)
    return inserted_nodes
end # function merge_trees!

"""
    check_node!(
        ref_node::T, leaves::Vector{T},
        leaf_ranks_reverse::Dict{String, Int64},
        xleft_dict::Dict{T, Tuple{T, Dict{Int64, T}}},
        xright_dict::Dict{T, Tuple{T, Dict{Int64, T}}},
        marked_nodes::Union{Dict{Int64, Bool, Nothing}}=nothing
        )::Union{Nothing, Tuple{T, T, T, Dict{Int64, T}, Dict{Int64, T}}} where T<:AbstractNode

--- INTERNAL ---
Helper function that handles major chunk of code that one_way_compatible and
merge_tree share
"""
function check_node!(ref_node::T, leaves::Vector{T},
                     leaf_ranks_reverse::Dict{String, Int64},
                     xleft_dict::Dict{T, Tuple{T, Dict{Int64, T}}},
                     xright_dict::Dict{T, Tuple{T, Dict{Int64, T}}},
                     marked_nodes::Union{Dict{Int64, Bool}, Nothing}=nothing
                     )::Union{Nothing, Tuple{T, T, T, T, T, Dict{Int64, T}, Dict{Int64, T}}} where T<:AbstractNode
    if ref_node.nchild != 0
        start = min_leaf_rank(leaf_ranks_reverse, ref_node)
        stop = max_leaf_rank(leaf_ranks_reverse, ref_node)
        start_node = leaves[start]
        stop_node = leaves[stop]
        xleft, left_path = xleft_dict[start_node]
        xright, right_path = xright_dict[stop_node]
        if intersect(values(left_path), values(right_path)) == []
            try marked_nodes[ref_node.num] = true catch; end # try
            return
        end # if
        !xleft.root && delete!(left_path, node_depth(xleft.mother))
        !xright.root && delete!(right_path, node_depth(xright.mother))
        depth_left = node_depth(xleft)
        depth_right = node_depth(xright)
        depth = depth_left >= depth_right
        depth ? p2 = right_path[depth_left] : p2 = left_path[depth_right]
        depth ? p1 = xleft : p1 = xright
        if p2 == p1
            r = p1
            if depth
                d = left_path[depth_left + 1]
                e = right_path[depth_left + 1]
            else
                d = left_path[depth_right + 1]
                e = right_path[depth_right + 1]
            end # if/else
        elseif p1.mother == p2.mother
            r = p1.mother
            depth ? d = p1 : d = p2
            depth ? e = p2 : e = p1
        else
            try marked_nodes[ref_node.num] = true catch; end # try
            return
        end # if/else
        if !(node_depth(d.mother) <= node_depth(r) && node_depth(e.mother) <= node_depth(r))
            try marked_nodes[ref_node.num] = true catch; return end # try
        else
            try marked_nodes[ref_node.num] = false
            catch
                isnothing(marked_nodes) && return (d, e, r, start_node, stop_node, left_path, right_path)
            end # try
            return
        end # if/else
    end # if
end # check_node

"""
    get_cluster_start_indeces(ref_nodes::T, tree::T)
        ::Dict{T, Int64} where T<:AbstractNode

--- INTERNAL ---
Helper function to obtain the cluster start indeces for a tree (tree), based
on the nodes of another tree (ref_nodes).
"""
function get_cluster_start_indeces(ref_nodes::Vector{T}, tree::T
                                  )::Dict{T, Int64} where T<:AbstractNode

    leaf_ranks = get_leaf_ranks(ref_nodes)
    leaves = Vector{Node}()
    cluster_start_indeces = Dict{Node, Int64}()
    nodes = post_order(tree)
    for node in nodes
        if node.nchild != 0
            for leaf in leaves
                path = split(node.binary, ",")
                if node_depth(node) <= node_depth(leaf) && path == split(leaf.binary, ",")[1:length(path)]
                    cluster_start_indeces[node] = leaf_ranks[leaf.name]
                    break
                end # if
            end # for
        else
            cluster_start_indeces[node] = leaf_ranks[node.name]
            push!(leaves, node)
        end # if / else
    end # for
    return cluster_start_indeces
end

"""
    function get_leaf_ranks(nodes::Vector{T})
        ::Dict{String, Int64} where T<:AbstractNode

--- INTERNAL ---
Enumerate the leaf nodes in a tree. Returns a dictionary of this mapping.
"""
function get_leaf_ranks(nodes::Vector{T}
                       )::Dict{String, Int64} where T<:AbstractNode

    leaf_ranks = Dict{String, Int64}()
    count = 1
    for node in nodes
        if node.nchild == 0
            leaf_ranks[node.name] = count
            count += 1
        end # if
    end # for
    return leaf_ranks
end # get_leaf_ranks

"""
    order_tree!(root::T, cluster_start_indeces::Dict{T, Int64}, leaves=Vector{T}())
        ::Vector{T} where T<:AbstractNode

--- INTERNAL ---
Helper function to order a tree based on cluster indeces and return the leaves
of the ordered tree
"""
function order_tree!(root::T, cluster_start_indeces::Dict{T, Int64},
                     leaves=Vector{T}())::Vector{T} where T<:AbstractNode

    sort!(root.children, by = child -> cluster_start_indeces[child],
          alg=InsertionSort)
    for child in root.children
        if child.nchild == 0
            push!(leaves, child)
        else
            order_tree!(child, cluster_start_indeces, leaves)
        end # if/else
    end # for
    set_binary!(root)
    return leaves
end # order_tree!

"""
    min_leaf_rank(leaf_ranks::Dict{String, Int64}, node::T)
        ::Int64 where T<:AbstractNode

---- INTERNAL ---
Recursive helper function to find the lowest ranked leaf descendant of a node
"""
function min_leaf_rank(leaf_ranks::Dict{String, Int64}, node::T
                      )::Int64 where T<:AbstractNode

    if node.nchild == 0
        return leaf_ranks[node.name]
    else
        possible_minima = Vector{Int64}()
        for child in node.children
            if child.nchild == 0
                push!(possible_minima, leaf_ranks[child.name])
            else
                push!(possible_minima, min_leaf_rank(leaf_ranks, child))
            end # if/else
        end # for
    end # if/else
    minimum(possible_minima)
end # min_leaf_rank

"""
    max_leaf_rank(leaf_ranks::Dict{String, Int64}, node::T)
        ::Int64 where T<:AbstractNode

--- INTERNAL ---
Recursive helper function to find the highest ranked leaf descendant of a node
"""
function max_leaf_rank(leaf_ranks::Dict{String, Int64}, node::T
                      )::Int64 where T<:AbstractNode

    if node.nchild == 0
        return leaf_ranks[node.name]
    else
        possible_maxima = Vector{Int64}()
        for child in node.children
            if child.nchild == 0
                push!(possible_maxima, leaf_ranks[child.name])
            else
                push!(possible_maxima, max_leaf_rank(leaf_ranks, child))
            end # if/else
        end # for
    end # if/else
    maximum(possible_maxima)
end # max_leaf_rank

"""
    x_left(node::T)::Tuple{T,Vector{T}} where T<:AbstractNode

--- INTERNAL ---
Helper function to find ancestor of a leaf that has said leaf as leftmost
descendant. Also returns the path from the leaf to the mother of that node.
"""
function x_left(node::T)::Tuple{T, Vector{T}} where T<:AbstractNode
    path = [node]
    while true
        if node.root
            return node, path
        else
            mother = node.mother
            if mother.children[1] != node
                push!(path, node.mother)
                return node, path
            end # if
        node = mother
        push!(path, node)
        end # if/else
    end # while
end # x_left

"""
    x_right(node::T)::Tuple{T,Vector{T}} where T<:AbstractNode

--- INTERNAL ---
Helper function to find ancestor of a leaf that has said leaf as rightmost
descendant. Also returns the path from the leaf to the mother of that node.
"""
function x_right(node::T)::Tuple{T, Vector{T}} where T<:AbstractNode
    path = [node]
    while true
        if node.root
            return node, path
        else
            mother = node.mother
            if mother.children[end] != node
                push!(path, node.mother)
                return node, path
            end # if
        node = mother
        push!(path, node)
        end # if/else
    end # while
end # x_right

"""
    depth_dicts(leaves::Vector{T})
        ::Tuple{Dict{T, Tuple{T, Dict{Int64, T}}}, Dict{T, Tuple{T, Dict{Int64, T}}}} where T<:AbstractNode

--- INTERNAL ---
Helper function for one_way_compatible and merge_trees!. Creates a dictionary
based on a vector of leaf nodes, and stores the depth of each node, as well as
the left and right path leading to it. Based on section 6.1 of the paper.
"""
function depth_dicts(leaves::Vector{T})::Tuple{Dict{T, Tuple{T, Dict{Int64, T}}}, Dict{T, Tuple{T, Dict{Int64, T}}}} where T<:AbstractNode
    xleft_dict = Dict{Node, Tuple{Node, Dict{Int64,Node}}}()
    xright_dict = Dict{Node, Tuple{Node, Dict{Int64, Node}}}()
    for leaf in leaves
        x, path = x_left(leaf)
        xleft_dict[leaf] = (x, Dict(node_depth(node) => node for node in path))
        x, path = x_right(leaf)
        xright_dict[leaf] = (x, Dict(node_depth(node) => node for node in path))
    end
    return xleft_dict, xright_dict
end


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the compressible ideal GLM-MHD equations
equations = IdealGlmMhdEquations2D(1.4)

initial_condition = initial_condition_weak_blast_wave

surface_flux = (flux_hindenlang_gassner, flux_nonconservative_powell)
volume_flux  = (flux_hindenlang_gassner, flux_nonconservative_powell)
polydeg = 4
basis = LobattoLegendreBasis(polydeg)
indicator_sc = IndicatorHennemannGassner(equations, basis,
                                         alpha_max=0.5,
                                         alpha_min=0.001,
                                         alpha_smooth=true,
                                         variable=density_pressure)
volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
                                                 volume_flux_dg=volume_flux,
                                                 volume_flux_fv=surface_flux)
solver = DGSEM(basis, surface_flux, volume_integral)

# Get the curved quad mesh from a mapping function
# Mapping as described in https://arxiv.org/abs/2012.12040
function mapping(xi, eta)
    y = 2.0 * eta + 1.0 / 6.0 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta))

    x = 2.0 * xi + 1.0 / 6.0 * (cos(0.5 * pi * xi) * cos(2 * pi * y))

    return SVector(x, y)
  end

cells_per_dimension = (16, 16)

mesh = StructuredMesh(cells_per_dimension, mapping, periodicity=true)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 0.4)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

cfl = 1.0
stepsize_callback = StepsizeCallback(cfl=cfl)

glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback,
                        stepsize_callback,
                        glm_speed_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

ext: jl


text:
# Code for Video 2.2:

fac(n) = (n <= 1)? 1 : n * fac(n - 1)

n = 20
@printf("The %d-th Factorial number is %d", n, fac(n))
println("")
print("The calculation took")
@time fac(n)

# The 20-th Factorial number is 2432902008176640000
# The calculation took  0.000006 seconds (149 allocations: 10.167 KB)

# C:\Users\CVO>julia --lisp
# ;  _
# ; |_ _ _ |_ _ |  . _ _
# ; | (-||||_(_)|__|_)|_)
# ;-------------------|------

# > (* 3 45)
# 135

# >

 n = 108 # 108
 n # 108
 :n # :n

 typeof(:n) # Symbol

 e = :(a + b) # :(a + b)

 typeof(e) # Expr

 a = 2 # 2
 b = 3 # 3
 eval(e) # 5

 fieldnames(e)
# 3-element Array{Symbol,1}:
#  :head
#  :args
#  :typ

 e.args
# 3-element Array{Any,1}:
#  :+
#  :a
#  :b

 c = quote
           a = 2
           b = 3
           a + b
       end
# quote  # none, line 2:
#     a = 2 # none, line 3:
#     b = 3 # none, line 4:
#     a + b
# end

 eval(c) # 5

 typeof(c) # Expr

 c = quote
           a = 2
           b = 3
           $a + $b
       end
# quote  # none, line 2:
#     a = 2 # none, line 3:
#     b = 3 # none, line 4:
#     2 + 3
# end

macro timeit(ex)
    quote
        local t0 = time()
        local val = $(esc(ex))
        local t1 = time()
        println("elapsed time in seconds: ")
        @printf "%.3f\n" t1 - t0
        val
    end
end

@timeit fac(20)
# elapsed time in seconds:
# 0.002
# 2432902008176640000

macro repeat(n, body)
    quote
        for i = 1:$(esc(n))
            $(esc(body))
        end
    end
end

@repeat(7, println("Hi Julia"))
# Hi Julia
# Hi Julia
# Hi Julia
# Hi Julia
# Hi Julia
# Hi Julia
# Hi Julia

macroexpand(:@repeat(7, println("Hi Julia")))
# quote  # none, line 3:
#     for #57#i = 1:7 # none, line 4:
#         println("Hi Julia")
#     end
# end

@repeat(3, for i in 1:5 println(i*i) end)
# 1
# 4
# 9
# 16
# 25
# 1
# 4
# 9
# 16
# 25
# 1
# 4
# 9
# 16
# 25

macroexpand(:@repeat(3, for i in 1:5 println(i*i) end))
# quote  # none, line 3:
#     for #29#i = 1:3 # none, line 4:
#         for i = 1:5 # none, line 1:
#             println(i * i)
#         end
#     end
# end

# @assert  108 = (2 + 105)
# ERROR: syntax: invalid assignment location "108"

# @assert  108 == (2 + 105)
# ERROR: AssertionError: 108 == 2 + 105

 using Base.Test

# @test_approx_eq 2 2.1
# ERROR: assertion failed: |2 - 2.1| <= 4.440892098500626e-12
#   2 = 2
#   2.1 = 2.1
#   difference = 0.10000000000000009 > 4.440892098500626e-12
#  in error at error.jl:22

# @test_approx_eq_eps 2 2.1 0.1
# ERROR: assertion failed: |2 - 2.1| <= 0.1
#   2 = 2
#   2.1 = 2.1
#   difference = 0.10000000000000009 > 0.1
#  in error at error.jl:22

@test_approx_eq_eps 2 2.1 0.15

@windows ? println("I am on Windows") : println("No Windows here")
# I am on Windows

@which (1 +2im) * (2 - 8im)
# *(z::Complex{T<:Real}, w::Complex{T<:Real}) at complex.jl:113

max_stars_repo_path: Section 2/2_2.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 2/2_2.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 2/2_2.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
using Test, CSV, Mmap, Dates, Tables, PooledArrays, CodecZlib, FilePathsBase, SentinelArrays, Parsers, WeakRefStrings, InlineStrings

const dir = joinpath(dirname(pathof(CSV)), "..", "test", "testfiles")

@eval macro $(:try)(ex)
    quote
        try $(esc(ex))
        catch
        end
    end
end

@testset "CSV" begin

@testset "CSV.File" begin

include("basics.jl")
include("testfiles.jl")
include("iteration.jl")

end # @testset "CSV.File"

include("write.jl")

@testset "PooledArrays" begin

    f = CSV.File(IOBuffer("X\nb\nc\na\nc"), pool=true)
    @test typeof(f.X) == PooledArrays.PooledArray{InlineString1,UInt32,1,Array{UInt32,1}}
    @test (length(f), length(f.names)) == (4, 1)
    @test f.X == ["b", "c", "a", "c"]
    @test f.X.refs[2] == f.X.refs[4]

    f = CSV.File(IOBuffer("X\nb\nc\na\nc"), pool=0.75)
    @test typeof(f.X) == PooledArrays.PooledArray{InlineString1,UInt32,1,Array{UInt32,1}}
    @test (length(f), length(f.names)) == (4, 1)
    @test f.X == ["b", "c", "a", "c"]
    @test f.X.refs[2] == f.X.refs[4]

    f = CSV.File(IOBuffer("X\nb\nc\n\nc"), pool=true, ignoreemptyrows=false)
    @test typeof(f.X) == PooledArray{Union{Missing, InlineString1},UInt32,1,Array{UInt32,1}}
    @test (length(f), length(f.names)) == (4, 1)
    @test f.X[3] === missing

    f = CSV.File(IOBuffer("X\nc\nc\n\nc\nc\nc\nc\nc\nc"), pool=0.25, ignoreemptyrows=false)
    @test typeof(f.X) == PooledArray{Union{Missing, InlineString1},UInt32,1,Array{UInt32,1}}
    @test (length(f), length(f.names)) == (9, 1)
    @test isequal(f.X, ["c", "c", missing, "c", "c", "c", "c", "c", "c"])

end

@testset "CSV.Rows" begin

    rows = CSV.Rows(IOBuffer("X\nb\nc\na\nc"))
    show(rows)
    row = first(rows)
    @test row.X == "b"
    @test row[1] == "b"

    # limit
    rows = collect(CSV.Rows(IOBuffer("X\nb\nc\na\nc"), limit=1))
    @test length(rows) == 1
    @test rows[1][1] == "b"
    rows = collect(CSV.Rows(IOBuffer("X\nb\nc\na\nc"), limit=0))
    @test length(rows) == 0

    # transpose
    rows = collect(CSV.Rows(IOBuffer("x,1\nx2,2\n"), transpose=true))
    @test length(rows) == 1
    @test rows[1].x == "1"
    @test rows[1].x2 == "2"

    # ignorerepeated
    rows = collect(CSV.Rows(IOBuffer("x   y   z\n   1   2   3"), ignorerepeated=true, delim=' '))
    @test length(rows) == 1
    @test length(rows[1]) == 3
    @test all(rows[1] .== ["1", "2", "3"])
    row = rows[1]
    @test CSV.Parsers.parse(Int, row, 1) == 1
    @test CSV.Parsers.parse(Int, row, :x) == 1
    @test CSV.detect(row, 2) == 2
    @test CSV.detect(row, :y) == 2

    # 448
    @test_throws ArgumentError CSV.Rows(IOBuffer("x\n1\n2\n3\n#4"), ignorerepeated=true)

    # 447
    rows = collect(CSV.Rows(IOBuffer("a,b,c\n1,2,3\n\n")))
    @test length(rows) == 1
    @test all(rows[1] .== ["1", "2", "3"])

    # not enough columns
    rows = collect(CSV.Rows(IOBuffer("x,y,z\n1\n2,3\n4,5,6")))
    @test length(rows) == 3
    @test all(isequal.(rows[1], ["1", missing, missing]))
    @test all(isequal.(rows[2], ["2", "3", missing]))
    @test all(rows[3] .== ["4", "5", "6"])

    # too many columns
    rows = collect(CSV.Rows(IOBuffer("x,y,z\n1,2,3\n4,5,6,7,8,")))
    @test length(rows) == 2
    @test isequal(values(rows[1]), ["1", "2", "3", missing, missing])
    @test isequal(values(rows[2]), ["4", "5", "6", "7", "8"])

    # fatal
    @test_throws CSV.Error collect(CSV.Rows(IOBuffer("x\n\"invalid quoted field")))

    # reusebuffer
    rows = collect(CSV.Rows(IOBuffer("x\n1\n2\n3")))
    @test length(rows) == 3
    @test rows[1][1] == "1"
    @test rows[2][1] == "2"
    @test rows[3][1] == "3"
    rows = collect(CSV.Rows(IOBuffer("x\n1\n2\n3"), reusebuffer=true))
    @test length(rows) == 3
    @test rows[1][1] == "3"
    @test rows[2][1] == "3"
    @test rows[3][1] == "3"

    for (i, row) in enumerate(CSV.Rows(IOBuffer("x\n1\n2\n3"), reusebuffer=true))
        @test all(row .== [string(i)])
    end

    # 903
    rows = collect(CSV.Rows(IOBuffer("int,float,date,datetime,bool,null,str,catg,int_float\n1,3.14,2019-01-01,2019-01-01T01:02:03,true,,hey,abc,2\n,,,,,,,,\n2,NaN,2019-01-02,2019-01-03T01:02:03,false,,there,abc,3.14\n"); types=[Int, Float64, Date, DateTime, Bool, Missing, String, String, Float64]))
    @test isequal(collect(rows[1]), [1, 3.14, Date(2019, 1, 1), DateTime(2019, 1, 1, 1, 2, 3), true, missing, "hey", "abc", 2.0])
    foreach(x -> @test(x === missing), rows[2])
    @test isequal(collect(rows[3]), [2, NaN, Date(2019, 1, 2), DateTime(2019, 1, 3, 1, 2, 3), false, missing, "there", "abc", 3.14])
end

@testset "CSV.detect" begin

@test CSV.detect("") == ""
@test CSV.detect("1") == 1
@test CSV.detect("1.1") == 1.1
@test CSV.detect("2015-01-01") == Date(2015)
@test CSV.detect("2015-01-01T03:04:05") == DateTime(2015, 1, 1, 3, 4, 5)
@test CSV.detect("03:04:05") == Time(3, 4, 5)
@test CSV.detect("true") === true
@test CSV.detect("false") === false
@test CSV.detect("abc") === "abc"

end

@testset "CSV.promote_types" begin

@test CSV.promote_types(Int64, Float64) == Float64
@test CSV.promote_types(Int64, Int64) == Int64
@test CSV.promote_types(CSV.NeedsTypeDetection, Float64) == Float64
@test CSV.promote_types(Float64, CSV.NeedsTypeDetection) == Float64
@test CSV.promote_types(Float64, Missing) == Float64
@test CSV.promote_types(Float64, String) === String

end

@testset "CSV.File with select/drop" begin

csv = """
a,b,c,d,e
1,2,3,4,5
6,7,8,9,10
"""

f = CSV.File(IOBuffer(csv), select=[1, 3, 5])
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), select=[:a, :c, :e])
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), select=["a", "c", "e"])
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), select=[true, false, true, false, true])
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), select=(i, nm) -> i in (1, 3, 5))
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), select=Int[])
@test length(f) == 2
@test length(f[1]) == 0

f = CSV.File(IOBuffer(csv), select=[1, 2, 3, 4, 5])
@test length(f) == 2
@test length(f[1]) == 5

f = CSV.File(IOBuffer(csv), drop=[2, 4])
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), drop=[:b, :d])
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), drop=["b", "d"])
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), drop=[false, true, false, true, false])
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), drop=(i, nm) -> i in (2, 4))
@test f.a == [1, 6]
@test all(f[1] .== [1, 3, 5])
@test length(f[1]) == 3

f = CSV.File(IOBuffer(csv), drop=Int[])
@test length(f) == 2
@test length(f[1]) == 5

f = CSV.File(IOBuffer(csv), drop=[1, 2, 3, 4, 5])
@test length(f) == 2
@test length(f[1]) == 0

end

@testset "CSV.Chunks" begin

chunks = CSV.Chunks(transcode(GzipDecompressor, Mmap.mmap(joinpath(dir, "randoms.csv.gz"))); ntasks=2)
@test length(chunks) == 2
state = iterate(chunks)
f, st = state
@test length(f) == 35086
state = iterate(chunks, st)
f, st = state
@test length(f) == 34914

end

function strs(x::Vector, e=nothing)
    len = sum(x -> sizeof(coalesce(x, "")), x)
    data = Vector{UInt8}(undef, len)
    poslens = Vector{PosLen}(undef, length(x))
    pos = 1
    anymissing = false
    for (i, s) in enumerate(x)
        if s !== missing
            slen = sizeof(s)
            bytes = codeunits(s)
            copyto!(data, pos, bytes, 1, slen)
            poslens[i] = PosLen(pos, slen, false, e !== nothing && e in bytes)
            pos += slen
        else
            anymissing = true
            poslens[i] = PosLen(pos, 0, true, false)
        end
    end
    return PosLenStringVector{anymissing ? Union{Missing, PosLenString} : PosLenString}(data, poslens, something(e, 0x00))
end

arrtype(x::ChainedVector{T, AT}) where {T, AT} = AT

@testset "CSV.chaincolumns!" begin

mv = MissingVector(1)
vec = [big(1)]
svec = convert(SentinelArray{BigInt}, [big(1), missing])
# MissingVector -> Vector{CustomType}
@test arrtype(CSV.chaincolumns!(ChainedVector([mv]), vec)) == typeof(svec)
# Vector{CustomType} -> MissingVector
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), mv)) == typeof(svec)
# MissingVector -> SentinelVector{CustomType}
@test arrtype(CSV.chaincolumns!(ChainedVector([mv]), svec)) == typeof(svec)
# SentinelVector{CustomType} -> MissingVector
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), mv)) == typeof(svec)
# Vector{CustomType} -> SentinelVector{CustomType}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), svec)) == typeof(svec)
# SentinelVector{CustomType} -> Vector{CustomType}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), vec)) == typeof(svec)
# Vector{CustomType} -> Vector{CustomType}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), vec)) == Vector{BigInt}
# SentinelVector{CustomType} -> SentinelVector{CustomType}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), svec)) == typeof(svec)

vec = strs(["hey", "there", "sailor"])
svec = strs(["hey", "ho", missing])
# MissingVector -> PosLenStringVector{PosLenString}
@test arrtype(CSV.chaincolumns!(ChainedVector([mv]), vec)) <: PosLenStringVector{Union{PosLenString, Missing}}
# PosLenStringVector{PosLenString} -> MissingVector
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), mv)) <: PosLenStringVector{Union{PosLenString, Missing}}
# MissingVector -> PosLenStringVector{Union{PosLenString, Missing}}
@test arrtype(CSV.chaincolumns!(ChainedVector([mv]), svec)) <: PosLenStringVector{Union{PosLenString, Missing}}
# PosLenStringVector{Union{PosLenString, Missing}} -> MissingVector
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), mv)) <: PosLenStringVector{Union{PosLenString, Missing}}
# PosLenStringVector{PosLenString} -> PosLenStringVector{Union{PosLenString, Missing}}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), svec)) <: PosLenStringVector{Union{PosLenString, Missing}}
# PosLenStringVector{Union{PosLenString, Missing}} -> PosLenStringVector{PosLenString}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), vec)) <: PosLenStringVector{Union{PosLenString, Missing}}

vec = [true, false, true]
svec = [true, false, missing]
# MissingVector -> Vector{Bool}
@test arrtype(CSV.chaincolumns!(ChainedVector([mv]), vec)) == Vector{Union{Bool, Missing}}
# Vector{Bool} -> MissingVector
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), mv)) == Vector{Union{Bool, Missing}}
# MissingVector -> Vector{Union{Bool, Missing}}
@test arrtype(CSV.chaincolumns!(ChainedVector([mv]), svec)) == Vector{Union{Bool, Missing}}
# Vector{Union{Bool, Missing}} -> MissingVector
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), mv)) == Vector{Union{Bool, Missing}}
# Vector{Bool} -> Vector{Union{Bool, Missing}}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), svec)) == Vector{Union{Bool, Missing}}
# Vector{Union{Bool, Missing}} -> Vector{Bool}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), vec)) == Vector{Union{Bool, Missing}}

vec = Int32[1, 2, 3]
svec = Union{Int32, Missing}[1, 2, missing]
# MissingVector -> Vector{Int32}
@test arrtype(CSV.chaincolumns!(ChainedVector([mv]), vec)) == Vector{Union{Int32, Missing}}
# Vector{Int32} -> MissingVector
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), mv)) == Vector{Union{Int32, Missing}}
# MissingVector -> Vector{Union{Int32, Missing}}
@test arrtype(CSV.chaincolumns!(ChainedVector([mv]), svec)) == Vector{Union{Int32, Missing}}
# Vector{Union{Int32, Missing}} -> MissingVector
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), mv)) == Vector{Union{Int32, Missing}}
# Vector{Int32} -> Vector{Union{Int32, Missing}}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec]), svec)) == Vector{Union{Int32, Missing}}
# Vector{Union{Int32, Missing}} -> Vector{Int32}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec]), vec)) == Vector{Union{Int32, Missing}}

vec1 = Int32[1, 2, 3]
vec2 = Int64[1, 2, 3]
svec1 = convert(SentinelArray{Int32}, [Int32(1), Int32(2), missing])
svec2 = convert(SentinelArray{Int64}, [Int64(1), Int64(2), missing])
# Vector{Int32} -> Vector{Int64}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec1]), vec2)) == Vector{Int64}
# Vector{Int64} -> Vector{Int32}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec2]), vec1)) == Vector{Int64}
# Vector{Int32} -> SentinelVector{Int64}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec1]), svec2)) == typeof(svec2)
# SentinelVector{Int32} -> Vector{Int64}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec1]), vec2)) == typeof(svec2)
# SentinelVector{Int32} -> SentinelVector{Int64}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec1]), svec2)) == typeof(svec2)
# SentinelVector{Int64} -> SentinelVector{Int32}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec2]), svec1)) == typeof(svec2)

vec1 = String3["hey", "ho"]
vec2 = String7["hey", "ho"]
svec1 = convert(SentinelArray{String3}, [String3("hey"), missing])
svec2 = convert(SentinelArray{String7}, [String7("hey"), missing])
# Vector{InlineString3} -> Vector{InlineString7}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec1]), vec2)) == Vector{String7}
# Vector{InlineString7} -> Vector{InlineString3}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec2]), vec1)) == Vector{String7}
# Vector{InlineString3} -> SentinelVector{InlineString7}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec1]), svec2)) == typeof(svec2)
# SentinelVector{InlineString3} -> Vector{InlineString7}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec1]), vec2)) == typeof(svec2)
# SentinelVector{InlineString3} -> SentinelVector{InlineString7}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec1]), svec2)) == typeof(svec2)
# SentinelVector{InlineString7} -> SentinelVector{InlineString3}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec2]), svec1)) == typeof(svec2)

vec1 = Int64[1, 2]
vec2 = [3.14, 2.28]
svec1 = convert(SentinelArray{Int64}, [Int64(1), missing])
svec2 = convert(SentinelArray{Float64}, [3.14, missing])
# Vector{Int64} -> Vector{Float64}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec1]), vec2)) == Vector{Float64}
# Vector{Float64} -> Vector{Int64}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec2]), vec1)) == Vector{Float64}
# Vector{Int64} -> SentinelVector{Float64}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec1]), svec2)) == typeof(svec2)
# SentinelVector{Int64} -> Vector{Float64}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec1]), vec2)) == typeof(svec2)
# SentinelVector{Int64} -> SentinelVector{Float64}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec1]), svec2)) == typeof(svec2)
# SentinelVector{Float64} -> SentinelVector{Int64}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec2]), svec1)) == typeof(svec2)

vec1 = String3["hey", "ho"]
vec2 = ["hey", "ho"]
svec1 = convert(SentinelArray{String3}, [String3("hey"), missing])
svec2 = convert(SentinelArray{String}, ["ho", missing])
# Vector{InlineString3} -> Vector{String}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec1]), vec2)) == Vector{String}
# Vector{String} -> Vector{InlineString3}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec2]), vec1)) == Vector{String}
# Vector{InlineString3} -> SentinelVector{String}
@test arrtype(CSV.chaincolumns!(ChainedVector([vec1]), svec2)) == typeof(svec2)
# SentinelVector{InlineString3} -> Vector{String}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec1]), vec2)) == typeof(svec2)
# SentinelVector{InlineString3} -> SentinelVector{String}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec1]), svec2)) == typeof(svec2)
# SentinelVector{String} -> SentinelVector{InlineString3}
@test arrtype(CSV.chaincolumns!(ChainedVector([svec2]), svec1)) == typeof(svec2)

vec1 = String3["hey", "ho"]
vec2 = ["hey", "ho"]
svec1 = convert(SentinelArray{String3}, [String3("hey"), missing])
svec2 = convert(SentinelArray{String}, ["ho", missing])
function pvecs(vec1, vec2, svec1, svec2)
    return PooledArray(vec1), PooledArray(vec2), PooledArray(svec1), PooledArray(svec2)
end
# MissingVector -> PooledVector{String7}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(mv, pvec1)) == typeof(psvec1)
# PooledVector{String} -> MissingVector
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(pvec2, mv)) == typeof(psvec2)
# MissingVector -> PooledVector{Union{String, Missing}}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(mv, psvec2)) == typeof(psvec2)
# PooledVector{Union{String7, Missing}} -> MissingVector
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(psvec1, mv)) == typeof(psvec1)
# PooledVector{String} -> Vector{String}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(pvec2, vec2)) == typeof(pvec2)
# PooledVector{String} -> SentinelVector{String7}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(pvec2, svec1)) == typeof(psvec2)
# PooledVector{Union{String, Missing}} -> Vector{String}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(psvec2, vec2)) == typeof(psvec2)
# PooledVector{Union{String7, Missing}} -> SentinelVector{String}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(psvec1, svec2)) == typeof(psvec2)
# PooledVector{String} -> PooledVector{Union{String, Missing}}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(pvec2, psvec2)) == typeof(psvec2)
# PooledVector{String} -> PooledVector{Union{String7, Missing}}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(pvec2, psvec1)) == typeof(psvec2)
# Vector{String} -> PooledVector{String7}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(vec2, pvec1)) == typeof(pvec2)
# SentinelVector{String} -> PooledVector{String7}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(svec2, pvec1)) == typeof(psvec2)
# Vector{String7} -> PooledVector{Union{String, Missing}}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(vec1, psvec2)) == typeof(psvec2)
# SentinelVector{String7} -> PooledVector{Union{String, Missing}}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(svec1, psvec2)) == typeof(psvec2)
# PooledVector{Union{String7, Missing}} -> PooledVector{String7}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(psvec1, pvec1)) == typeof(psvec1)
# PooledVector{Union{String, Missing}} -> PooledVector{String7}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(psvec2, pvec1)) == typeof(psvec2)

# issue 925; ambiguity in `CSV._promote(::Type{PooledVector{T,R,RA}}, x::PooledVector{T,R,RA})`
# PooledVector{String3} -> PooledVector{String}
pvec1, pvec2, psvec1, psvec2 = pvecs(vec1, vec2, svec1, svec2)
@test typeof(CSV.chaincolumns!(pvec1, pvec2)) == typeof(pvec2)

end

@testset "CSV.File vector inputs" begin

data = [
    "a,b,c\n1,2,3\n4,5,6\n",
    "a,b,c\n7,8,9\n10,11,12\n",
    "a,b,c\n13,14,15\n16,17,18",
]

f = CSV.File(map(IOBuffer, data))
@test length(f) == 6
@test f.names == [:a, :b, :c]
@test f.a == [1, 4, 7, 10, 13, 16]

data = [
    "a,b\nbill,x\njane,y\n",
    "a,b\ntomm,z\ntimm,\n",
    "a,b\njoee,z\njerr,g\n",
]

f = CSV.File(map(IOBuffer, data))
@test length(f) == 6
@test isequal(f.b, ["x", "y", "z", missing, "z", "g"])

data = [
    "a,b,c\n1,2,3\n4,5,6\n",
    "a,b,c\n7.14,8,9\n10,11,12\n",
    "a,b,c\n13,14,15\n16,17,18",
]

f = CSV.File(map(IOBuffer, data))
@test eltype(f.a) == Float64

data = [
    "a,b,c\n1,2,3\n4,5,6\n",
    "a2,b,c\n7,8,9\n10,11,12\n",
    "a,b,c\n13,14,15\n16,17,18",
]

f = CSV.File(map(IOBuffer, data))
@test isequal(f.a, [1, 4, missing, missing, 13, 16])

f = CSV.File(map(IOBuffer, data); source=:source)
@test eltype(f.source) == String
@test f.source isa PooledArray

f = CSV.File(map(IOBuffer, data); source="source")
@test eltype(f.source) == String
@test f.source isa PooledArray

f = CSV.File(map(IOBuffer, data); source="source"=>[1,2,3])
@test eltype(f.source) <: Integer
@test f.source isa PooledArray

f = CSV.File(map(IOBuffer, data); source=:source=>["1", "2", "3"])
@test eltype(f.source) == String
@test f.source isa PooledArray

end

end


text:
using Rows
using Base.Test

A = [1 2; 3 4; 5 6]

@test [[1 2], [3 4], [5 6]] == rows(A)
@test [1 2] == rows(A)[1]
@test [[1 2], [3 4]] == rows(A)[1:2]
@test [[1 2], [3 4], [5 6]] == rows(A)[1:end]
@test [[1 2], [3 4], [5 6]] == [row for row in rows(A)]

@test [[1; 3; 5], [2; 4; 6]] == columns(A)
@test [1; 3; 5] == columns(A)[1]
@test [[1; 3; 5], [2; 4; 6]] == columns(A)[1:end]
@test [[1; 3; 5], [2; 4; 6]] == [col for col in columns(A)]

B = ['1' '2'; '3' '4'; '5' '6'; '7' '8']
@test [['1' '2'], ['3' '4'], ['5' '6'], ['7' '8']] == rows(B)

C = ["1" "2"; "3" "4"; "5" "6"; "7" "8"]
@test [["1" "2"], ["3" "4"], ["5" "6"], ["7" "8"]] == rows(C)


text:
using Conda, Pkg, Test, VersionParsing

exe = Sys.iswindows() ? ".exe" : ""

Conda.update()

if Conda.USE_MINIFORGE
    default_channel = "conda-forge"
    alt_channel = "defaults"
else
    default_channel = "defaults"
    alt_channel = "conda-forge"
end

env = :test_conda_jl
rm(Conda.prefix(env); force=true, recursive=true)

@testset "Install Python package" begin
    Conda.add("python", env)
    pythonpath = joinpath(Conda.python_dir(env), "python" * exe)
    @test isfile(pythonpath)

    cmd = Conda._set_conda_env(`$pythonpath -c "import six"`, env)
    @test_throws Exception run(cmd)
    Conda.add("six", env)
    run(cmd)
end

@testset "Install executable package" begin
    @test Conda.exists("curl", env)
    Conda.add("curl", env)

    curlvers = Conda.version("curl",env)
    @test curlvers >= v"1.0"
    @test Conda.exists("curl==$curlvers", env)

    curl_path = joinpath(Conda.bin_dir(env), "curl" * exe)
    @test isfile(curl_path)

    @test "curl" in Conda.search("cu*", env)

    Conda.rm("curl", env)
    @test !isfile(curl_path)
end

pythonpath = joinpath(Conda.PYTHONDIR, "python" * exe)
@test isfile(pythonpath)
pyversion = read(`$pythonpath -c "import sys; print(sys.version)"`, String)
@test pyversion[1:1] == Conda.MINICONDA_VERSION

condarc_path = joinpath(homedir(), ".condarc")
last_modified = mtime(condarc_path)

Conda.add_channel("foo", env)
@test Conda.channels(env) == ["foo", default_channel]

# Testing that calling the function twice do not fail
Conda.add_channel("foo", env)

# Validate that only the Conda.jl RC file was modified
@test occursin("foo", read(Conda.conda_rc(env), String))
@test !isfile(condarc_path) || !occursin("foo", read(condarc_path, String))

Conda.rm_channel("foo", env)
@test mtime(condarc_path) == last_modified

@test Conda.channels(env) == [default_channel]

# Add a package from a specific channel
Conda.add("zlib", env; channel=alt_channel)

@testset "Batch install and uninstall" begin
    Conda.add(["affine", "ansi2html"], env)
    installed = Conda._installed_packages(env)
    @test "affine" ∈ installed
    @test "ansi2html" ∈ installed

    Conda.rm(["affine", "ansi2html"], env)
    installed = Conda._installed_packages(env)
    @test "affine" ∉ installed
    @test "ansi2html" ∉ installed
end

# Run conda clean
Conda.clean(; debug=true)

@testset "Exporting and creating environments" begin
    new_env = :test_conda_jl_2
    Conda.add("curl", env)
    Conda.export_list("conda-pkg.txt", env)

    # Create a new environment
    rm(Conda.prefix(new_env); force=true, recursive=true)
    Conda.import_list(
        IOBuffer(read("conda-pkg.txt")), new_env; channels=["foo", alt_channel, default_channel]
    )

    # Ensure that our new environment has our channels and package installed.
    @test Conda.channels(new_env) == ["foo", alt_channel, default_channel]
    installed = Conda._installed_packages(new_env)
    @test "curl" ∈ installed
    rm("conda-pkg.txt")
end

@testset "Conda.pip_interop" begin
    Conda.pip_interop(false, env)
    @test_throws Exception Conda.check_pip_interop(env)
    @test !Conda.pip_interop(env)

    Conda.pip_interop(true, env)
    @test Conda.pip_interop(env)
    @test Conda.check_pip_interop(env)
end

@testset "Conda.pip" begin
    Conda.pip("install", "affine", env)
    @test "affine" ∈ Conda._installed_packages(env)
    Conda.pip("uninstall -y", "affine", env)
    @test "affine" ∉ Conda._installed_packages(env)

    Conda.pip("install", ["affine", "ansi2html"], env)
    installed = Conda._installed_packages(env)
    @test "affine" ∈ installed
    @test "ansi2html" ∈ installed

    Conda.pip("uninstall -y", ["affine", "ansi2html"], env)
    installed = Conda._installed_packages(env)
    @test "affine" ∉ installed
    @test "ansi2html" ∉ installed
end

@testset "build" begin
    condadir = abspath(first(DEPOT_PATH), "conda")
    depsfile = joinpath(@__DIR__, "..", "deps", "deps.jl")

    function preserve_build(body)
        condadir_bak = condadir * "-backup"
        depsfile_bak = depsfile * "-backup"

        # Backup contents that may be changed by a build.
        ispath(condadir) && mv(condadir, condadir_bak)
        ispath(depsfile) && mv(depsfile, depsfile_bak)

        try
            body()
        finally
            # Restore content
            ispath(condadir_bak) && mv(condadir_bak, condadir; force=true)
            ispath(depsfile_bak) && mv(depsfile_bak, depsfile; force=true)
        end
    end

    if Sys.ARCH in [:x86, :i686]
        CONDA_JL_USE_MINIFORGE_DEFAULT = "false"
    else
        CONDA_JL_USE_MINIFORGE_DEFAULT = "true"
    end

    @testset "defaults" begin
        preserve_build() do
            # In order to test the defaults no depsfiles must be present
            @test !isfile(depsfile)
            @test !isfile(joinpath(condadir, "deps.jl"))

            withenv("CONDA_JL_VERSION" => nothing, "CONDA_JL_HOME" => nothing, "CONDA_JL_USE_MINIFORGE" => nothing) do
                Pkg.build("Conda")
                @test read(depsfile, String) == """
                    const ROOTENV = "$(escape_string(joinpath(condadir, "3")))"
                    const MINICONDA_VERSION = "3"
                    const USE_MINIFORGE = $(CONDA_JL_USE_MINIFORGE_DEFAULT)
                    """
            end
        end
    end

    @testset "miniforge" begin
        preserve_build() do
            # In order to test the defaults no depsfiles must be present
            @test !isfile(depsfile)
            @test !isfile(joinpath(condadir, "deps.jl"))

            withenv("CONDA_JL_VERSION" => nothing, "CONDA_JL_HOME" => nothing, "CONDA_JL_USE_MINIFORGE" => "1") do
                Pkg.build("Conda")
                @test read(depsfile, String) == """
                    const ROOTENV = "$(escape_string(joinpath(condadir, "3")))"
                    const MINICONDA_VERSION = "3"
                    const USE_MINIFORGE = true
                    """
            end
        end
        preserve_build() do
            # In order to test the defaults no depsfiles must be present
            @test !isfile(depsfile)
            @test !isfile(joinpath(condadir, "deps.jl"))

            withenv("CONDA_JL_VERSION" => nothing, "CONDA_JL_HOME" => nothing, "CONDA_JL_USE_MINIFORGE" => "0") do
                Pkg.build("Conda")
                @test read(depsfile, String) == """
                    const ROOTENV = "$(escape_string(joinpath(condadir, "3")))"
                    const MINICONDA_VERSION = "3"
                    const USE_MINIFORGE = false
                    """
            end
        end
    end

    @testset "custom home" begin
        preserve_build() do
            mktempdir() do dir
                withenv("CONDA_JL_VERSION" => "3", "CONDA_JL_HOME" => dir, "CONDA_JL_USE_MINIFORGE" => nothing) do
                    Pkg.build("Conda")
                    @test read(depsfile, String) == """
                        const ROOTENV = "$(escape_string(dir))"
                        const MINICONDA_VERSION = "3"
                        const USE_MINIFORGE = $(CONDA_JL_USE_MINIFORGE_DEFAULT)
                        """
                end
            end
        end
    end

    @testset "version mismatch" begin
        preserve_build() do
            # Mismatch in written file
            write(depsfile, """
                const ROOTENV = "$(escape_string(joinpath(condadir, "3")))"
                const MINICONDA_VERSION = "2"
                const USE_MINIFORGE = $(CONDA_JL_USE_MINIFORGE_DEFAULT)
                """)

            withenv("CONDA_JL_VERSION" => nothing, "CONDA_JL_HOME" => nothing, "CONDA_JL_USE_MINIFORGE" => nothing) do
                Pkg.build("Conda")
                @test read(depsfile, String) == """
                    const ROOTENV = "$(escape_string(joinpath(condadir, "2")))"
                    const MINICONDA_VERSION = "2"
                    const USE_MINIFORGE = $(CONDA_JL_USE_MINIFORGE_DEFAULT)
                    """
            end

            # ROOTENV should be replaced since CONDA_JL_HOME wasn't explicitly set
            withenv("CONDA_JL_VERSION" => "3", "CONDA_JL_HOME" => nothing, "CONDA_JL_USE_MINIFORGE" => nothing) do
                Pkg.build("Conda")
                @test read(depsfile, String) == """
                    const ROOTENV = "$(escape_string(joinpath(condadir, "3")))"
                    const MINICONDA_VERSION = "3"
                    const USE_MINIFORGE = $(CONDA_JL_USE_MINIFORGE_DEFAULT)
                    """
            end
        end
    end
end


text:
using GeometricIntegrators.Utils

function rel_energy_err(sol::SolutionSDE{AT,TT,WT,1}) where {AT,TT,WT}
    en_ref  = 0.5 * L2norm(sol.q[begin])
    en_last = 0.5 * L2norm(sol.q[end])

    return abs((en_last-en_ref)/en_ref)
end

function rel_energy_err(sol::SolutionSDE{AT,TT,WT,2}) where {AT,TT,WT}
    en_ref  = 0.5 .* map(q -> L2norm(q), sol.q[begin,:])
    en_last = 0.5 .* map(q -> L2norm(q), sol.q[end,:])

    return maximum(abs.( (en_last .- en_ref) ./ en_ref ))
end

function rel_energy_err(sol::SolutionPSDE{AT,TT,WT,1}) where {AT,TT,WT}
    en_ref  = 0.5 * ( L2norm(sol.q[begin]) + L2norm(sol.p[begin]) )
    en_last = 0.5 * ( L2norm(sol.q[end])   + L2norm(sol.p[end])   )

    return abs((en_last-en_ref)/en_ref)
end

function rel_energy_err(sol::SolutionPSDE{AT,TT,WT,2}) where {AT,TT,WT}
    en_ref  = 0.5 .* map(q -> L2norm(q), sol.q[begin,:]) .+
              0.5 .* map(p -> L2norm(p), sol.p[begin,:])
    en_last = 0.5 .* map(q -> L2norm(q), sol.q[end,:]) .+
              0.5 .* map(p -> L2norm(p), sol.p[end,:])

    return maximum(abs.( (en_last .- en_ref) ./ en_ref ))
end


text:
using neutrinos
using Test

@testset "neutrinos.jl" begin

    #----- Basic tests for computing the matter mixing angle

    # Test basic vacuum angles. Start with with theta = 0, s.t. tan(2*theta)
    @test 0.0 == neutrinos.getMatterMixingAngle(0.0,0.0);

    ## If the coupling strength is zero, the denominator is 1.0. At 2theta = pi/4.0, the numerator is also 1 such that the result is theta_m = pi/8, which correctly matches the vacuum angle.
    @test isapprox(neutrinos.getMatterMixingAngle(pi/8.0,0.0), pi/8.0; atol=1.0e-5);

    # There is an analytical result nearly as simple as the above, but with coupling strength equal to 1.0. In this case, the denominator is 1 - sqrt(2.0).
    @test isapprox(neutrinos.getMatterMixingAngle(pi/8.0,1.0), 0.5*atan(1.0/(1.0-sqrt(2.0))); atol=1.0e-5);

    ## This version exploits the sqrt(3.0) produced by cos(pi/6.0) (which is just cos(2theta) for theta = pi/12).

    #Really not sure where this one came from...

    @test isapprox(neutrinos.getMatterMixingAngle(pi/12.0,sqrt(3.0)), 0.5*atan(-0.5); atol=1.0e-5);

    ## Problem 12.6 in Guidry's stars book presents an alternative formulation for the mixing angle based on sin that can be used to make the comparsion here.
    twoTheta = 1.0;
    chi = 0.5;
    thetaM_fromSin = 0.5*asin(sin(twoTheta)/(sqrt(1.0-2.0*chi*cos(twoTheta)+chi*chi)));
    @test isapprox(thetaM_fromSin,neutrinos.getMatterMixingAngle(twoTheta/2.0,chi); atol=1.0e-10);

    #----- Basic tests for computing the coupling strength

    # Test zero case. Yeah, this means that the neutrinos have no energy and that there are no electrons in the universe! :-P
    @test 0.0 == neutrinos.getCouplingStrength(0.0,0.0);

    # Test the case where either the neutrino energy is zero or the electron number density is zero. Both cases result in no coupling strength.
    @test 0.0 == neutrinos.getCouplingStrength(0.0,5.0);
    @test 0.0 == neutrinos.getCouplingStrength(5.0,0.0);

    # If the energy is equal to deltaMSquared/(2sqrt(2)*fermiConst), then the coupling strength is equal to the electron number density.
    @test isapprox( neutrinos.getCouplingStrength((neutrinos.deltaMSquared/(2.0*sqrt(2.0)*neutrinos.fermiConst)),1.0),1.0; atol=1.0e-5);

    # If the energy and the electron density are both one, the coupling strength is just 2sqrt(2)fermiConst/deltaMSquared.
    @test isapprox( neutrinos.getCouplingStrength(1.0,1.0), (2.0*sqrt(2.0)*neutrinos.fermiConst/neutrinos.deltaMSquared); atol=1.0e-5);

    # I made this one up using the defaults and r/R = 0.0 for the number density. It uses chi = L/l_m directly.
    l_m = sqrt(2.0)*pi/(neutrinos.fermiConst*1.475e26);
    L = 4.0*pi*neutrinos.energy/neutrinos.deltaMSquared;
    chi = L/l_m;
    println("test: $L, $l_m, $chi");
    @test isapprox(chi,neutrinos.getCouplingStrength(neutrinos.energy,1.475e26); atol=1.0e-5);

    #----- Tests for the scaling factor f(theta,chi)
    # For chi = theta = 0.0, f = 1.
    @test isapprox(1.0,neutrinos.getScalingFactor(0.0,0.0);atol=1.0e-16)
    # For chi = 1.0, theta = 5.0pi/6.0, f = 1.0
    @test isapprox(1.0,neutrinos.getScalingFactor(5.0*pi/6.0,1.0);atol=1.0e-8)
    # For chi = 1.0, theta = 2.0pi/6.0, f = sqrt(3.0)
    @test isapprox(sqrt(3.0),neutrinos.getScalingFactor(2.0*pi/6.0,1.0);atol=1.0e-8)
    # For chi = 1.0, theta = 4.0pi/6.0, f = sqrt(3.0)
    @test isapprox(sqrt(3.0),neutrinos.getScalingFactor(4.0*pi/6.0,1.0);atol=1.0e-8)

    #----- Basic tests for the electron number density
    # This simple test checks behavior when R/R_sun = 0.0, which sets the exponent to zero and the result is the coefficient, 1.475E26.
    @test isapprox(1.475E26, neutrinos.getElectronNumberDensity(0.0),atol=1.0e-5);
    # Similar to the above, but in the case where the radius is equal to the solar radius and the fraction is 1.0, then the result is 1.475E26*exp(-10.54).
    @test isapprox(1.475E26*exp(-10.54), neutrinos.getElectronNumberDensity(1.0), atol=1.0e-5);

    #----- Basic tests for the electron neutrino probability
    # For r/R = theta = 0.0, P = 1.0.
    @test 1.0 == neutrinos.getElectronNeutrinoProbability(0.0,0.0);
    #@test 0.5 == #neutrinos.getElectronNeutrinoProbability();
    # This is the example from Guidry's stars book (ADD DEETS!)
end


text:
"""

Package `StructuredArrays` implements arrays whose elements have values that
only depend on the indices.  Example of such arrays are uniform arrays with
constant values, structured arrays with boolean values indicating whether an
entry is significant.

"""
module StructuredArrays

export
    AbstractStructuredArray,
    AbstractStructuredMatrix,
    AbstractStructuredVector,
    AbstractUniformArray,
    AbstractUniformMatrix,
    AbstractUniformVector,
    MutableUniformArray,
    MutableUniformMatrix,
    MutableUniformVector,
    StructuredArray,
    StructuredMatrix,
    StructuredVector,
    UniformArray,
    UniformMatrix,
    UniformVector

using ArrayTools
import Base: @propagate_inbounds

abstract type AbstractStructuredArray{T,N,S<:IndexStyle} <:
    AbstractArray{T,N} end

abstract type AbstractUniformArray{T,N} <:
    AbstractStructuredArray{T,N,IndexLinear} end

"""
    UniformArray(val, siz) -> A

yields an array `A` which behaves as an immutable array of size `siz` whose
values are all equal to `val`.  The storage requirement is `O(1)` instead of
`O(prod(siz))` for a usual array.  The array dimensions may be specified as
multiple arguments.

Uniform arrays implement conventional linear indexing: `A[i]` yields `val`
for all linear indices `i` in the range `1:length(A)`.

A statement like `A[i] = val` is not implemented as uniform arrays are
considered as immutable.  Call `MutableUniformArray(val,siz)` to create a
uniform array whose element value can be changed.

""" UniformArray

struct UniformArray{T,N} <: AbstractUniformArray{T,N}
    len::Int
    siz::Dims{N}
    val::T
    UniformArray{T,N}(val, siz::Dims{N}) where {T,N} =
        new{T,N}(checksize(siz), siz, val)
end

"""
    MutableUniformArray(val, siz) -> A

yields an array `A` which behaves as a mutable array of size `siz` whose
values are all `val`.  The storage requirement is `O(1)` instead of
`O(prod(siz))` for a usual array.  The array dimensions may be specified as
multiple arguments.

Uniform arrays implement conventional linear indexing: `A[i]` yields `val`
for all linear indices `i` in the range `1:length(A)`.

A statement like `A[i] = val` is allowed but changes the value of all the
elements of `A`.  Call `UniformArray(val,siz)` to create an immutable uniform
array whose element value cannot be changed.

""" MutableUniformArray

mutable struct MutableUniformArray{T,N} <: AbstractUniformArray{T,N}
    len::Int
    siz::Dims{N}
    val::T
    MutableUniformArray{T,N}(val, siz::Dims{N}) where {T,N} =
        new{T,N}(checksize(siz), siz, val)
end

"""
    StructuredArray([S = IndexCartesian,] fnc, siz) -> A

yields an array `A` which behaves as an array of size `siz` whose values are a
given function, here `fnc`, of its indices: `A[i]` is computed as `fnc(i)`.
The storage requirement is `O(1)` instead of `O(prod(siz))` for a usual array.
The array dimensions may be specified as multiple arguments.

The optional argument `S` may be used to specifiy another index style than the
default `IndexCartesian`, for instance `IndexLinear`.  If specified, `S` may be
a sub-type of `IndexStyle` or an instance of such a sub-type.  If `S` is
`IndexCartesian` (the default), the function `fnc` will be called with `N`
integer arguments, `N` being the number of dimensions.  If `S` is
`IndexCartesian`, the function `fnc` will be called with a single integer
argument.

For instance, the structure of a lower triangular matrix of size `m×n` would be
given by:

    StructuredArray((i,j) -> (i ≥ j), m, n)

but with a constant small storage requirement whatever the size of the matrix.

Although the callable object `fnc` may not be a *pure function*, its return
type shall be stable and structured arrays are considered as immutable in the
sense that a statement like `A[i] = val` is not implemented.  The type of the
elements of structured array is guessed by applying `fnc` to the unit index.
The element type, say `T`, may also be explicitely specified:

    StructuredArray{T}([S = IndexCartesian,] fnc, siz)

""" StructuredArray

struct StructuredArray{T,N,S,F} <: AbstractStructuredArray{T,N,S}
    len::Int
    siz::Dims{N}
    fnc::F
    StructuredArray{T,N,S,F}(fnc, siz::Dims{N}) where {T,N,S,F} =
        new{T,N,S,F}(checksize(siz), siz, fnc)
end

const AbstractStructuredVector{T,S} = AbstractStructuredArray{T,1,S}
const StructuredVector{T,S,F} = StructuredArray{T,1,S,F}
const AbstractUniformVector{T} = AbstractUniformArray{T,1}
const MutableUniformVector{T} = MutableUniformArray{T,1}
const UniformVector{T} = UniformArray{T,1}

const AbstractStructuredMatrix{T,S} = AbstractStructuredArray{T,2,S}
const StructuredMatrix{T,S,F} = StructuredArray{T,2,S,F}
const AbstractUniformMatrix{T} = AbstractUniformArray{T,2}
const MutableUniformMatrix{T} = MutableUniformArray{T,2}
const UniformMatrix{T} = UniformArray{T,2}

# Specialize base abstract array methods for StructuredArray, UniformArray, and
# MutableUniformArray and provide basic constructors to convert trailing
# arguments to dimensions.
for cls in (:StructuredArray, :UniformArray, :MutableUniformArray)
    @eval begin
        Base.length(A::$cls) = getfield(A, :len)
        Base.size(A::$cls) = getfield(A, :siz)
        Base.size(A::$cls{T,N}, i::Integer) where {T,N} =
            (i < 1 ? bad_dimension_index() : i ≤ N ? @inbounds(size(A)[i]) : 1)
        Base.axes1(A::$cls{T,0}) where {T} = Base.OneTo(1)
        Base.axes1(A::$cls{T,N}) where {T,N} = Base.OneTo(size(A)[1])
        Base.axes(A::$cls) = map(Base.OneTo, size(A))
        Base.axes(A::$cls, i::Integer) = Base.OneTo(size(A, i))
        Base.has_offset_axes(::$cls) = false

        $cls(arg1, dims::Integer...) = $cls(arg1, dims)
        $cls{T}(arg1, dims::Integer...) where {T} = $cls{T}(arg1, dims)
        $cls{T,N}(arg1, dims::Integer...) where {T,N} = $cls{T,N}(arg1, dims)
    end
end
for cls in (:UniformArray, :MutableUniformArray)
    @eval begin
        $cls(val::T, siz::NTuple{N,Integer}) where {T,N} =
            $cls{T,N}(val, to_size(siz))
        $cls{T}(val, siz::NTuple{N,Integer}) where {T,N} =
            $cls{T,N}(val, to_size(siz))
        $cls{T,N}(val, siz::NTuple{N,Integer}) where {T,N} =
            $cls{T,N}(val, to_size(siz))
    end
end

StructuredArray(fnc, siz::NTuple{N,Integer}) where {N} =
    StructuredArray(IndexCartesian, fnc, to_size(siz))
StructuredArray{T}(fnc, siz::NTuple{N,Integer}) where {T,N} =
    StructuredArray{T}(IndexCartesian, fnc, to_size(siz))
StructuredArray{T,N}(fnc, siz::NTuple{N,Integer}) where {T,N} =
    StructuredArray{T}(fnc, siz)

# All constructors for StructuredArray are based on the 2 first ones (only
# depending on whether parameter T is provided or not) so the other constructors
# just make sure and arguments have correct type.

function StructuredArray(::Type{S},
                         fnc::F,
                         siz::Dims{N}) where {N,S<:IndexStyle,F}
    T = guess_eltype(fnc, S, Val(N))
    StructuredArray{T,N,S,F}(fnc, siz)
end

function StructuredArray{T}(::Type{S},
                            fnc::F,
                            siz::Dims{N}) where {T,N,S<:IndexStyle,F}
    StructuredArray{T,N,S,F}(fnc, siz)
end

guess_eltype(fnc, ::Type{IndexLinear}, ::Val{N}) where {N} =
    typeof(fnc(one(Int)))

guess_eltype(fnc, ::Type{IndexCartesian}, ::Val{N}) where {N} =
    typeof(fnc(ntuple(i -> one(Int), Val(N))...))

# Index style specified and size specified as a tuple or by trailing arguments.

function StructuredArray(::Union{S,Type{S}}, fnc,
                         siz::Integer...) where {S<:IndexStyle}
    StructuredArray(S, fnc, siz)
end
function StructuredArray(::Union{S,Type{S}}, fnc,
                         siz::NTuple{N,Integer}) where {N,S<:IndexStyle}
    StructuredArray(S, fnc, to_size(siz))
end

function StructuredArray{T}(::Union{S,Type{S}}, fnc,
                            siz::Integer...) where {T,S<:IndexStyle}
    StructuredArray{T}(S, fnc, siz)
end
function StructuredArray{T}(::Union{S,Type{S}}, fnc,
                            siz::NTuple{N,Integer}) where {T,N,S<:IndexStyle}
    StructuredArray{T}(S, fnc, to_size(siz))
end

function StructuredArray{T,N}(::Union{S,Type{S}}, fnc,
                              siz::Integer...) where {T,N,S<:IndexStyle}
    StructuredArray{T,N}(S, fnc, siz) # keep the N to check
end
function StructuredArray{T,N}(::Union{S,Type{S}}, fnc,
                              siz::NTuple{N,Integer}) where {T,N,S<:IndexStyle}
    StructuredArray{T}(S, fnc, to_size(siz))
end

function StructuredArray{T,N,S}(::Union{S,Type{S}}, fnc,
                                siz::Integer...) where {T,N,S<:IndexStyle}
    StructuredArray{T,N}(S, fnc, siz) # keep the N to check
end
function StructuredArray{T,N,S}(::Union{S,Type{S}}, fnc,
                                siz::NTuple{N,Integer}) where {T,N,S<:IndexStyle}
    StructuredArray{T}(S, fnc, to_size(siz))
end

function StructuredArray{T,N,S}(fnc,
                                siz::Integer...) where {T,N,S<:IndexStyle}
    StructuredArray{T,N}(S, fnc, siz) # keep the N to check
end
function StructuredArray{T,N,S}(fnc,
                                siz::NTuple{N,Integer}) where {T,N,S<:IndexStyle}
    StructuredArray{T}(S, fnc, to_size(siz))
end

# See comments near `getindex` in `abstractarray.jl` for explanations about how
# `getindex` and `setindex!` methods are expected to be specialized depending
# on the indexing style.

Base.IndexStyle(::Type{<:AbstractStructuredArray{T,N,S}}) where {T,N,S} = S()

@inline function Base.getindex(A::StructuredArray{T,N,IndexLinear},
                               i::Int) :: T where {T,N}
    @boundscheck checkbounds(A, i)
    A.fnc(i)
end

@inline function Base.getindex(A::StructuredArray{T,N,IndexCartesian},
                               inds::Vararg{Int, N}) :: T where {T,N}
    @boundscheck checkbounds(A, inds...)
    A.fnc(inds...)
end

@inline function Base.getindex(A::X,
                               i::Int) where {X<:Union{<:UniformArray,
                                                       <:MutableUniformArray}}
    @boundscheck checkbounds(A, i)
    A.val
end

@inline function Base.setindex!(A::MutableUniformArray, x, i::Int)
    @boundscheck checkbounds(A, i)
    (i == length(A) == 1) || not_all_elements()
    A.val = x
    A
end

@inline function Base.setindex!(A::MutableUniformArray, x,
                                i::AbstractUnitRange{<:Integer})
    (first(i) == 1 && last(i) == length(A)) || not_all_elements()
    A.val = x
    A
end

@inline function Base.setindex!(A::MutableUniformArray, x, ::Colon)
    A.val = x
    A
end

"""
    checksize(dims) -> len

yields the number of elements of an array of size `dims` throwing an error if
any dimension is invalid.

"""
function checksize(dims::Dims{N}) where {N}
    len = 1
    @inbounds for i in 1:N
        (dim = dims[i]) ≥ 0 || bad_dimension_length()
        len *= dim
    end
    return len
end

@noinline not_all_elements() =
    error("all elements must be set at the same time")

@noinline bad_dimension_index() = error("out of range dimension index")

@noinline bad_dimension_length() = error("invalid dimension length")

end # module


text:
"""
Structure used to track the context - the expression being modified and the path to the current node
"""
mutable struct GrammarContext
	originalExpr::RuleNode    	# original expression being modified
	nodeLocation::Vector{Int}   	# path to he current node in the expression, a sequence of child indices for each parent
end

GrammarContext(originalExpr::RuleNode) = GrammarContext(originalExpr, [])

"""
add parent to the context
"""
function addparent!(context::GrammarContext, parent::Int)
	push!(context.nodeLocation, parent)
end

"""
	Copies the given context and insert the parent in the node location
"""
function copy_and_insert(old_context::GrammarContext, parent::Int)
	new_context = GrammarContext(old_context.originalExpr, deepcopy(old_context.nodeLocation))
	push!(new_context.nodeLocation, parent)
	new_context
end


text:
"""
$(TYPEDEF)

Default prior on all dimensions of a variable node in the factor graph.  `Prior` is
not recommended when non-Euclidean dimensions are used in variables.
"""
struct Prior{T <: SamplableBelief} <: AbstractPrior
  Z::T
end
Prior(::UniformScaling) = Prior(Normal())

getManifold(pr::Prior) = TranslationGroup(getDimension(pr.Z))
# getSample(cf::CalcFactor{<:Prior}) = rand(cf.factor.Z)

# basic default
(s::CalcFactor{<:Prior})(z, x1) = z .- x1

## packed types are still developed by hand.  Future versions would likely use a @packable macro to write Protobuf safe versions of factors

"""
$(TYPEDEF)

Serialization type for Prior.
"""
Base.@kwdef mutable struct PackedPrior <: AbstractPackedFactor
  Z::PackedSamplableBelief
end
function convert(::Type{PackedPrior}, d::Prior)
  PackedPrior(convert(PackedSamplableBelief, d.Z))
end
function convert(::Type{Prior}, d::PackedPrior)
  Prior(convert(SamplableBelief, d.Z))
end

#


text:
module TDMSReader

import BitOperations:bget
import DataStructures:OrderedDict

include("types.jl")

const _example_tdms = (@__DIR__) * "\\..\\test\\example_files\\reference_file.tdms"
const _example_incremental = [(@__DIR__) * "\\..\\test\\example_files\\incremental_test_$i.tdms" for i = 1:6]
const _example_DAQmx = (@__DIR__) * "\\..\\test\\example_files\\DAQmx example.tdms"

export readtdms

readtdms() = readtdms(_example_tdms)
function readtdms(fn::AbstractString)
    s = open(fn)
    f = File()
    objdict = ObjDict()
    while !eof(s)
        temp = readleadin(s)
        isnothing(temp) && break
        toc, nextsegmentoffset, rawdataoffset = temp
        toc.kTocNewObjList && empty!(objdict.current)
        meta_start = position(s)
        toc.kTocMetaData && readmetadata!(f, objdict, s)
        raw_len = nextsegmentoffset - rawdataoffset
        toc.kTocRawData && raw_len > 0 && readrawdata!(objdict, raw_len, seek(s, meta_start + rawdataoffset))
        seek(s, meta_start + nextsegmentoffset)
    end
    close(s)
    f
end

function readleadin(s::IO)
    if ntoh(read(s, UInt32)) != 0x5444_536D
        @warn "head mismatch"
        return nothing
    end
    toc = ToC(ltoh(read(s, UInt32)))
    toc.kTocBigEndian && throw(ErrorException("Big Endian files not supported"))
    read(s, UInt32) == 4713 || throw(ErrorException("File not recongnized as TDMS formatted file"))
    (toc = toc, nextsegmentoffset = read(s, UInt64), rawdataoffset = read(s, UInt64),)
end

function readmetadata!(f::File, objdict::ObjDict, s::IO)
    n = read(s, UInt32)
    for i = 1:n
        readobj!(f, objdict, s)
    end
end

read(s::IO, ::Type{T}) where {T} = Base.read(s, T)
read(s::IO, ::Type{String}) where {T} = String(Base.read(s, read(s, UInt32)))
function readobj!(f::File, objdict::ObjDict, s::IO)
    objpath = read(s, String)
    println(objpath)
    # @info "Read @ position $(hexstring(position(s)))"
    rawdata = read(s, UInt32)
    hasrawdata = hasnewchunk = false
    if rawdata == 0xFF_FF_FF_FF # No Raw Data
    elseif iszero(rawdata) # Keep Chunk layout
        haskey(objdict.full, objpath) || throw(ErrorException("Previous Segment Missing"))
        hasrawdata = true
    elseif rawdata in (0x00_00_12_69, 0x00_00_13_69)
        readDAQmx(s::IO, rawdata)
        throw(ErrorException("Not Implemented"))
    else
        hasrawdata = hasnewchunk = true
        T = tdsTypes[read(s, UInt32)]
        T <: TDMSUnimplementedType && throw(ErrorException("TDMS Data Type of $T is not supported"))
        T == String && throw(ErrorException("Need Functionality to Read String as raw data"))
        read(s, UInt32) == 1 || throw(ErrorException("TDMS Array Dimension is not 1"))
        n = read(s, UInt64)
        println("rawtype:", T, ", len:", n)
    end

    if objpath == "/"
        hasrawdata && throw(ErrorException("TDMS root should not have raw data"))
        readprop!(f.props, s)
    else
        m = match(r"\/'(.+?)'(?:\/'(.+)')?", objpath)
        isnothing(m) && throw(ErrorException("Object Path $objpath is malformed"))
        group, channel = m.captures
        g = get!(Returns(Group()), f, group)
        if isnothing(channel) # Is a Group
            hasrawdata && throw(ErrorException("TDMS Group should not have raw data"))
            readprop!(g.props, s)
        else # Is  a Channel
            chan = get!(g, channel) do
                hasrawdata ? Channel{T}() : Channel{Nothing}()
            end
            readprop!(chan.props, s)
        end
    end
    if hasrawdata
        if hasnewchunk
            objdict.full[objpath] = Chunk{T}(chan.data, n)
        end
        objdict.current[objpath] = objdict.full[objpath]
    end
end

function readprop!(props, s::IO)
    for i = 1:read(s, UInt32)
        propname = read(s, String)
        T = tdsTypes[read(s, UInt32)]
        props[propname] = read(s, T)
    end
end

function readrawdata!(objdict::ObjDict, totalbytes::Integer, s::IO)
    current = values(objdict.current)
    nbytes = sum(current) do val
        sizeof(eltype(val.data)) * val.nsamples
    end
    batch::Int = fld(totalbytes, nbytes)
    nbytes * batch == totalbytes || throw("layout error")
    # Accelerate via preallocation
    datas = foreach(current) do val
        resize!(val.data, length(val.data) + val.nsamples * batch)
    end
    for i in batch-1:-1:0
        foreach(current) do val
            n = Int(val.nsamples)
            ind = length(val.data) .+ (1 - n:0)
            @views read!(s, val.data[ind .- i * n])
        end
    end
end

function seekalign(s::IO)
    x = position(s)
    mask = typeof(x)(0b11)
    seek(s, ifelse(x & mask > 0, x  & ~mask + 4, x))
end

function readDAQmx(s::IO, id)
    @info "Read DAQmc Raw Data @ $(hexstring(id))"
    T = tdsTypes[read(s, UInt32)]
    @info "Data type $T"
    read(s, UInt32) == 1 || throw(ErrorException("TDMS Array Dimension is not 1"))
    chunksize = read(s, UInt64)
    @info "Chunk size = $chunksize"

    @info "Read vector of format change scalers"
    vectorsize = read(s, UInt32)
    V = tdsTypes[read(s, UInt32)]
    rawbufferindex = read(s, UInt32)
    rawbyteoffset = read(s, UInt32)
    sampleformatbitmap = read(s, UInt32)
    scaleid = read(s, UInt32)
    @info "Vector Size = $vectorsize"
    @info "DAQmx data type = $V"
    @info "Raw Buffer Index = $rawbufferindex"
    @info "Raw byte offset = $rawbyteoffset"
    @info "Sample Format Bitmatp = $sampleformatbitmap"
    @info "Scale ID = $scaleid"

    for i in 1:n

    end

end

    function hexstring(x::Integer)
    "0x$(lpad(string(x, base=16), 8, '0'))"
end

function readseginfo(fn::AbstractString)
    # STILL COULD USE SOME WORK HERE TO MAKE IT CLEAN AND NEAT
    s = open(fn)
    f = File()
    objdict = ObjDict()
    tdmsSeg = OrderedDict()
    lead_size = 28 # lead in is 28 bytes
    segCnt = 0
    while !eof(s)
        startPos = position(s)
        temp = readleadin(s)
        isnothing(temp) && break
        toc, nextsegmentoffset, rawdataoffset = temp
        if toc.kTocNewObjList
            empty!(objdict.current)
        end

        nobj = read(s, UInt32)
        seek(s, startPos + lead_size)
        if toc.kTocMetaData
            readmetadata!(f, objdict, s)
        end
        if toc.kTocRawData
            # readrawdata!(objdict, nextsegmentoffset-rawdataoffset, s)
        end
        nextsegmentPos = Int64(startPos + nextsegmentoffset + lead_size)
        dataPos = Int64(startPos + lead_size + rawdataoffset)
        rawdatasize = nextsegmentoffset - rawdataoffset
        tdmsTmp = SegInfo(startPos, toc, nextsegmentPos, dataPos, nobj, rawdatasize)
        segCnt += 1
        tdmsSeg[string("seg", segCnt)] = tdmsTmp
        seek(s, nextsegmentPos)
    end
    close(s)
    return f, objdict, tdmsSeg
end

end # module


text: fib(n) = ([1 1 ; 1 0]^n)[1,2]

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode