text:
####################################################################
# KD Tree
####################################################################
# The KDTree type
immutable KDTree{T <: AbstractFloat}
    data::Matrix{T} # dim x n_p array with floats
    hyper_recs::Vector{HyperRectangle{T}} # Each hyper rectangle bounds its children
    indices::Vector{Int} # Translation between tree index to point indes or adsa
    split_vals::Vector{T} # what values we split the tree at for a given internal node
    split_dims::Vector{Int} # what dimension we split the tree for a given internal node
    n_d::Int # Number of dimension in tree
    last_size::Int # Number of points in last node
    leafsize::Int # Number of points in each node except the last node
    n_leafs::Int # Number of leaf nodes
    n_internal::Int # Number of internal nodes
    cross_node::Int # Index of first node in last row
    first_leaf_row::Int # The row which contains the first leaf node
    offset::Int # How many nodes in first_leaf_row before the first leaf node
    data_reordered::Bool # If the data has been reordered for optimizing memory access.
end

function show(io::IO, tree::KDTree)
    print(io, string("KDTree from ", size(tree.data, 2),
                     " points in ", size(tree.data, 1), " dimensions."))
end

# Helper functions to get node numbers and points
getleft(idx::Int) = 2idx
getright(idx::Int) = 2idx + 1
getparent(idx::Int) = div(idx, 2)
isleaf(tree::KDTree, idx::Int) = idx > tree.n_internal

# Returns the number of points in a node.
# All nodes have leafsize points in them except possibly the last node.
function n_ps(n_leafs::Int, n_internal::Int, leafsize::Int,
              last_size::Int, idx::Int)
    if idx != n_leafs + n_internal
        return leafsize
    else
        return last_size
    end
end

# Returns the index for the first point for a given leaf node.
function point_index(cross_node::Int, offset::Int, last_size:: Int,
                     leafsize::Int, n_internal::Int, idx::Int)
    if idx >= cross_node
        return (idx - cross_node) * leafsize + 1
    else
        return ((offset + idx - 1 - n_internal) * leafsize
                  + last_size + 1)
    end
end

# Returns the range of indices for any node, internal or terminal.
# For an internal node the range of indices are of all the points
# that that subtree contains.
# This range might be split into two subranges.
# The first two return values are the indices for the first range and the
# second two are the indices for the second range.
function node_indices(tree, index)
    if isleaf(tree, index)
        p_index = point_index(tree.cross_node, tree.offset, tree.last_size,
                              tree.leafsize, tree.n_internal, index)
        n_p =  n_ps(tree.n_leafs, tree.n_internal, tree.leafsize,
                    tree.last_size, index)
        l1 = p_index
        r1 = p_index + n_p -1
        return l1, r1, 0, -1
    else
        first_leaf = tree.n_internal + 1
        left, right = index, index
        # Find most left node
        while (left < first_leaf) left = getleft(left) end

        # Find most right node
        while (right < first_leaf) right = getright(right) end

        if left < right
            return node_indices(tree, left)[1], node_indices(tree, right)[2], 0, -1
        else # here we need to use the subranges
            l1 = node_indices(tree, first_leaf)[1]
            r1 = node_indices(tree, right)[2]
            l2 = node_indices(tree, left)[1]
            r2 = node_indices(tree, tree.n_leafs + tree.n_internal)[2]
            return l1, r1, l2, r2
        end
    end
end

# Constructor for KDTree
function KDTree{T <: AbstractFloat}(data::Matrix{T};
                                    leafsize::Int = 10,
                                    reorder::Bool = true)

    if size(data, 2) == 0
        error("Need at least 1 point")
    end

    if size(data, 1) == 0
        error("Need at least 1 dimensional points")
    end

    n_d, n_p = size(data)

    n_leaf =  ceil(Integer, n_p / leafsize)
    n_internal = n_leaf - 1

    # Row whith first leaf node
    l = floor(Integer, log2(n_leaf))

    # Node index of first leaf node in last row
    cross_node = @compat Int(2^(l+1))

    # This only happens when n_p / leafsize is a power of 2?
    if cross_node >= n_internal + n_leaf
        cross_node = div(cross_node, 2)
    end

    # Number of nodes into row l where leaf nodes start
    offset = 2(n_leaf - 2^l) - 1

    # Last node might not have leafsize points in it.
    last_size = n_p % leafsize
    if last_size == 0
        last_size = leafsize
    end

    indices = collect(1:n_p)

    if reorder
        indices_reordered = Array(Int, n_p)
        data_reordered = Array(T, n_d, n_p)
    else
        # Dummy variables here
        indices_reordered = Array(Int, 0)
        data_reordered = Array(T, n_d, 0)
    end

    split_vals = Array(T, n_internal)
    split_dims = Array(Int, n_internal)
    hyper_recs = Array(HyperRectangle{T}, n_internal + n_leaf)

    # Create first bounding hyper rectangle
    maxes = Array(T, n_d)
    mins = Array(T, n_d)
    for j in 1:n_d
        dim_max = typemin(T)
        dim_min = typemax(T)
        for k in 1:n_p
            dim_max = max(data[j, k], dim_max)
            dim_min = min(data[j, k], dim_min)
        end
        maxes[j] = dim_max
        mins[j] = dim_min
    end

    hyper_recs[1] = HyperRectangle{T}(mins, maxes)

    low = 1
    high = n_p

    # Call the recursive KDTree builder
    build_KDTree(1, data, data_reordered, split_vals, split_dims, hyper_recs,
                 indices, indices_reordered, leafsize, low, high, cross_node,
                 offset, last_size, n_internal, reorder)

    if reorder
        KDTree(data_reordered, hyper_recs, indices_reordered, split_vals, split_dims, n_d, last_size,
              leafsize, n_leaf, n_internal, cross_node, l ,offset, reorder)
    else
        KDTree(data, hyper_recs, indices, split_vals, split_dims, n_d, last_size,
              leafsize, n_leaf, n_internal, cross_node, l ,offset, reorder)
    end

end

# Recursive function to build the tree.
# Calculates what dimension has the maximum spread,
# and how many points to send to each side.
# Sorts the indices at the split_
# Splits the hyper cubes and calls recursively
# with the new cubes and node indices.
# TODO: The number of arguments are growing ridiculous.
function build_KDTree{T <: AbstractFloat}(index::Int,
                                          data::Matrix{T},
                                          data_reordered::Matrix{T},
                                          split_vals::Vector{T},
                                          split_dims::Vector{Int},
                                          hyper_recs::Vector{HyperRectangle{T}},
                                          indices::Vector{Int},
                                          indices_reorder::Vector{Int},
                                          leafsize::Int,
                                          low::Int,
                                          high::Int,
                                          cross_node::Int,
                                          offset::Int,
                                          last_size::Int,
                                          n_internal::Int,
                                          reorder::Bool)

    n_p = high - low + 1 # Points left in this subtree

    if n_p <= leafsize

        if reorder
            # Here we reorder the data points so that points contained
            # in nodes with an index close to each other are also themselves
            # close in memory.
            p_index = point_index(cross_node, offset, last_size,
                                  leafsize, n_internal, index)
            n_p_node = n_ps(n_internal + 1, n_internal, leafsize,
                            last_size, index)

            for i in p_index:(p_index + n_p_node - 1)
                idx = indices[i]
                @devec data_reordered[:, i] = data[:, idx]

                # Saves the inverse n
                indices_reorder[i] = idx
            end
        end
        return
    end

    # The number of leafs left, ceil to count a partially filled node as 1.
    n_leafs = ceil(Integer, n_p / leafsize)

    # Rows left in the sub tree
    k = floor(Integer, log2(n_leafs))

    # Number of leftover nodes needed
    rest = n_leafs - 2^k

    # If the last leaf node will be on the right side of the tree we
    # send points so that left tree will be perfectly filled,
    # else we do the opposite.
    if k == 0
        mid_idx = low
    elseif n_p <= 2 * leafsize
        mid_idx = leafsize + low
    elseif  rest > 2^(k-1) # Last node over the "half line" in the row
        mid_idx = 2^k * leafsize + low
    elseif rest == 0 # Perfectly filling the last row
        mid_idx = 2^(k-1)* leafsize + low
    else
        mid_idx = n_p - 2^(k-1) * leafsize + low
    end

    n_d = size(data, 1)
    split_dim = 1
    max_spread = zero(T)
    hyper_rec = hyper_recs[index]
    # Find dimension and and spread where the spread is maximal
    for d in 1:n_d
        spread = hyper_rec.maxes[d] - hyper_rec.mins[d]
        if spread > max_spread
            max_spread = spread
            split_dim = d
        end
    end

    split_dims[index] = split_dim

    # sorts the points in the maximum spread dimension s.t
    # data[split_dim, a]) < data[split_dim, b]) for all a <= mid_idx, b > mid_idx
    _select!(indices, mid_idx, low, high, data, split_dim)

    split_val = data[split_dim, indices[mid_idx]]
    split_vals[index] = split_val

    # Create the hyper rectangles for the children
    hyper_rec_1, hyper_rec_2 = split_hyper_rec(hyper_recs[index], split_dim, split_val)
    hyper_recs[getleft(index)] = hyper_rec_1
    hyper_recs[getright(index)] = hyper_rec_2

    # Recursive call left and right sub tree
    build_KDTree(getleft(index), data, data_reordered,
                 split_vals, split_dims, hyper_recs,
                 indices, indices_reorder, leafsize, low, mid_idx - 1,
                 cross_node, offset, last_size,
                 n_internal, reorder)

    build_KDTree(getright(index), data, data_reordered,
                 split_vals, split_dims, hyper_recs,
                 indices, indices_reorder, leafsize, mid_idx, high,
                 cross_node, offset, last_size,
                 n_internal, reorder)
end

####################################################################
# Distances
####################################################################
# Reduced euclidian distances
@inline function euclidean_distance_red{T <: AbstractFloat}(point_1::AbstractVector{T},
                                                            point_2::AbstractVector{T})
    dist = zero(T)
    @simd for i in eachindex(point_1)
        @inbounds dist += abs2(point_1[i] - point_2[i])
    end
    return dist
end

@inline function euclidean_distance_red{T <: AbstractFloat}(data::Matrix{T},
                                                            idx::Int,
                                                            point::AbstractVector{T})
    dist = zero(T)
    @simd for i in eachindex(point)
        @inbounds dist += abs2(data[i, idx] - point[i])
    end
    return dist
end

@inline function euclidean_distance_red{T <: AbstractFloat}(tree::KDTree{T},
                                                            idx::Int,
                                                            tree2::KDTree{T},
                                                            idx2::Int)
    dist = 0.0
    for i = 1:tree.n_d
        @inbounds dist += abs2(tree.data[i, idx] - tree2.data[i, idx2])
    end
    return dist
end

####################################################################
# Query functions
####################################################################
# Finds the k nearest neighbour to a given point in space..
# The full_rec_dim argument requires some explanation.
# There is a tradeoff between improving the culling in the tree and
# the extra time that improvement take compared to just checking some
# extra nodes. If full_rec_dim > the dimensions in the tree the actual
# minimal distance between the rectangle and the point will be used for
# culling. Else, only the distance between the point and the rectangle
# for the dimension the rectangle was split at is used. As the dimension
# is increasing it is more and more important to use the actual distance.
# The tipping point seems to be right now for around a dimension of 6
# which is the value that will be used default.
function knn{T <: AbstractFloat}(tree::KDTree, point::Vector{T}, k::Int, full_rec_dim::Int = 6)

    # Check that k is not greater than points in tree
    if k > size(tree.data, 2) || k <= 0
        error("k > number of points in tree or <= 0")
    end

    # Check consistent dimension between tree and point
    if size(point,1) != size(tree.data, 1)
        error(string("Wrong dimension of input point, points in the tree",
                     " have dimension ", size(tree.data, 1), " you",
                     " gave a point with dimension ", size(point,1), "."))
    end

    # Initiate heaps to store indices and distances
    best_idxs = [-1 for _ in 1:k]
    best_dists = [typemax(T) for _ in 1:k]

    if tree.n_d < full_rec_dim
        _knn_small(tree, point, k, best_idxs, best_dists, 1)
    else
        init_min = get_min_distance(tree.hyper_recs[1], point)
        _knn(tree, point, k, best_idxs, best_dists, 1, init_min)
    end

    # Sqrt here because distances are stored in reduced format.
    @inbounds @simd for i in eachindex(best_dists)
        best_dists[i] = sqrt(best_dists[i])
    end

    # Translate indices back to the original if we reoredered the data
    if tree.data_reordered
        for j in 1:k
            @inbounds best_idxs[j] = tree.indices[best_idxs[j]]
        end
    end

    # Sort both heaps according to the distance.
    heap_sort_inplace!(best_dists, best_idxs)

    return best_idxs, best_dists
end

# Convert
knn{P <: Real}(tree::KDTree, point::Vector{P}, k::Int) =
  knn(tree, float(point), k)

function _knn{T <: AbstractFloat}(tree::KDTree{T},
                                  point::Vector{T},
                                  k::Int,
                                  best_idxs ::Vector{Int},
                                  best_dists::Vector{T},
                                  index::Int,
                                  min_dist::T)

    data_reordered = tree.data_reordered
    indices = tree.indices

    # If leaf, brute force through the points in the node and
    # if the distance is smaller add both the distance and index
    # to their respective heaps.
    if isleaf(tree, index)
        p_index = point_index(tree.cross_node, tree.offset, tree.last_size,
                              tree.leafsize, tree.n_internal, index)
        n_p =  n_ps(tree.n_leafs, tree.n_internal, tree.leafsize,
                    tree.last_size, index)
        @inbounds for z in p_index:p_index + n_p - 1
            idx = data_reordered ? z : indices[z]
            dist_d = euclidean_distance_red(tree.data, idx, point)
            if dist_d <= best_dists[1]
                best_dists[1] = dist_d
                best_idxs[1] = idx
                percolate_down!(best_dists, best_idxs, dist_d, idx)
            end
        end
        return
    end

    # Find what subtree is closest to the point
    if point[tree.split_dims[index]] < tree.split_vals[index]
        close_node = getleft(index)
        far_node = getright(index)
    else
        far_node = getleft(index)
        close_node = getright(index)
    end

    # Call close subtree first to improve culling
    _knn(tree, point, k, best_idxs, best_dists, close_node, min_dist)

    rec = tree.hyper_recs[close_node]
    dim = tree.split_dims[index]
    min_d_dim = get_min_dim(rec, point, dim)

    far_min = min_dist - min_d_dim + abs2(point[tree.split_dims[index]] - tree.split_vals[index])

    # Only call the sub tree further away if it is close enough
    if far_min < best_dists[1]
            _knn(tree, point, k, best_idxs, best_dists, far_node, far_min)
    end
    return
end

function _knn_small{T <: AbstractFloat}(tree::KDTree{T},
                                        point::Vector{T},
                                        k::Int,
                                        best_idxs ::Vector{Int},
                                        best_dists::Vector{T},
                                        index::Int)

    data_reordered = tree.data_reordered
    indices = tree.indices

    # If leaf, brute force through the points in the node and
    # if the distance is smaller add both the distance and index
    # to their respective heaps.
    if isleaf(tree, index)
        p_index = point_index(tree.cross_node, tree.offset, tree.last_size,
                              tree.leafsize, tree.n_internal, index)
        n_p =  n_ps(tree.n_leafs, tree.n_internal, tree.leafsize,
                    tree.last_size, index)
        @inbounds for z in p_index:p_index + n_p - 1
            idx = data_reordered ? z : indices[z]
            dist_d = euclidean_distance_red(tree.data, idx, point)
            if dist_d <= best_dists[1]
                best_dists[1] = dist_d
                best_idxs[1] = idx
                percolate_down!(best_dists, best_idxs, dist_d, idx)
            end
        end
        return
    end

    # Find what subtree is closest to the point
    if point[tree.split_dims[index]] < tree.split_vals[index]
        close_node = getleft(index)
        far_node = getright(index)
    else
        far_node = getleft(index)
        close_node = getright(index)
    end

    # Call close subtree first to improve culling
    _knn_small(tree, point, k, best_idxs, best_dists, close_node)

    # Only call the sub tree further away if it is close enough
    if abs2(point[tree.split_dims[index]] - tree.split_vals[index]) < best_dists[1]
            _knn_small(tree, point, k, best_idxs, best_dists, far_node)
    end
    return
end

# Returns the list of indices closer than a given point at a given
# radius.
function inball{T <: AbstractFloat}(tree::KDTree{T},
                                              p::Vector{T},
                                              radius::T,
                                              sort::Bool = false)

    # Check dimension consistency
    if size(p, 1) != size(tree.data, 1)
        error(string("Wrong dimension of input point, points in the tree",
                     " have dimension ", size(tree.data, 1), " you",
                     " gave a point with dimension ", size(p,1), "."))
    end

    index = 1
    idx_in_ball = Int[]

    init_min, init_max = get_min_max_distance(tree.hyper_recs[1], p)
    _in_ball(tree, index, p, abs2(radius) , idx_in_ball, init_min, init_max)

    if tree.data_reordered
        for i in 1:length(idx_in_ball)
            idx_in_ball[i] = tree.indices[idx_in_ball[i]]
        end
    end

    if sort
        sort!(idx_in_ball)
    end

    return idx_in_ball
end

inball{P <: Real, R <: Real}(tree::KDTree, p::Vector{P}, r::R) =
  inball(tree, float(p), float(r))

# Explicitly check the distance between leaf node and point while traversing
function _in_ball{T <: AbstractFloat}(tree::KDTree{T},
                                            index::Int,
                                            p::Vector{T},
                                            r::T,
                                            idx_in_ball::Vector{Int},
                                            min_dist::T,
                                            max_dist::T)

    if min_dist > r # Hyper sphere is outside hyper rectangle, skip the whole sub tree
        return
    end

    if max_dist < r
        addall(tree, index, idx_in_ball)
        return
    end

    if isleaf(tree, index)
        p_index = point_index(tree.cross_node, tree.offset, tree.last_size,
                              tree.leafsize, tree.n_internal, index)
        n_p =  n_ps(tree.n_leafs, tree.n_internal, tree.leafsize,
                    tree.last_size, index)
        for z in p_index:p_index + n_p - 1
            idx = tree.data_reordered ? z : tree.indices[z]
            dist_d = euclidean_distance_red(tree.data, idx, p)
            if dist_d < r
                push!(idx_in_ball, idx)
            end
        end
        return
    end

    d =  tree.split_dims[index]

    # We are querying 3 times so makes no sense to do this optimization unless n_d > 3.
    if tree.n_d > 3
        # # Remove contribution from this rectangle
        min_d_dim, max_d_dim = get_min_max_dim(tree.hyper_recs[index], p, d)
        max_dist -= max_d_dim
        min_dist -= min_d_dim

        # Add contribution of left node
        min_d_dim, max_d_dim = get_min_max_dim(tree.hyper_recs[getleft(index)], p, d)
        _in_ball(tree, getleft(index), p, r, idx_in_ball, min_dist + min_d_dim, max_dist + max_d_dim)

        # Add contribution of right node
        min_d_dim, max_d_dim = get_min_max_dim(tree.hyper_recs[getright(index)], p, d)
        _in_ball(tree, getright(index), p, r, idx_in_ball, min_dist + min_d_dim, max_dist + max_d_dim)
    else
        min_dist, max_dist = get_min_max_distance(tree.hyper_recs[index], p)
        _in_ball(tree, getleft(index), p, r, idx_in_ball, min_dist, max_dist)
        _in_ball(tree, getright(index), p, r, idx_in_ball, min_dist, max_dist)
    end

end

# Adds everything in this subtree since we have determined
# that the hyper rectangle completely encloses the hyper sphere
function addall(tree::KDTree, index::Int, idx_in_ball::Vector{Int})
    l1, r1, l2, r2 =  node_indices(tree, index)
    for r in (l1:r1, l2:r2)
        for z in r
            idx = tree.data_reordered ? z : tree.indices[z]
            push!(idx_in_ball, idx)
        end
    end
    return
end

# Returns the sorted list of indices for all points in the tree inside a
# hypersphere of a given point with a given radius.
function inball{T <: AbstractFloat}(tree::KDTree{T},
                                    tree2::KDTree{T},
                                    r::T,
                                    sort::Bool = false)

    if tree.data_reordered || tree2.data_reordered
        error(string("Trees with reoredered data is not currently supported.",
              " Please create your trees using for example KDTree(data, 10, true)"))
    end

    # Check tree dimensions consistency
    if size(tree.data, 1) != size(tree2.data, 1)
        error(string("Wrong dimension of input tree, points in the tree",
                     " have dimension ", size(tree.data, 1), " you",
                     " gave a tree with dimension ", size(tree2.data, 1), "."))
    end

    index = 1
    idx_in_ball =  Array((Array{Int64, 1}), size(tree.data, 2))
    for i in 1:length(idx_in_ball)
        idx_in_ball[i] = Int[]
    end

    init_min, init_max = get_min_max_distance(tree.hyper_recs[1], tree2.hyper_recs[1])

    # Square r since we work in reduced distance
    _in_ball(tree, index, tree2, index, abs2(r) , idx_in_ball, init_min, init_max)

    if sort
        for i in 1:length(idx_in_ball)
            sort!(idx_in_ball[i])
        end
    end

    return idx_in_ball
end

inball{P <: Real, R <: Real}(tree::KDTree, p::Vector{P}, r::R) =
  inball(tree, float(p), float(r))

# Explicitly check the distance between leaf node and point while traversing
function _in_ball{T <: AbstractFloat}(tree::KDTree{T},
                                      index::Int,
                                      tree2::KDTree{T},
                                      index2::Int,
                                      r::T,
                                      idx_in_ball::Vector{Vector{Int}},
                                      min_dist::T,
                                      max_dist::T)

    # This is currently ineffective for large d since we recalculate the distance for all dimensions
    # We should instead just calculate the change in distance for the current split_dim
    min_dist, max_dist = get_min_max_distance(tree.hyper_recs[index], tree2.hyper_recs[index2])
    if min_dist > r # Hyper sphere is outside hyper rectangle, skip the whole sub tree
        return

    elseif max_dist < r
        addall(tree, index, tree2, index2, idx_in_ball)

    elseif isleaf(tree, index)
        if isleaf(tree2, index2)
            l1, r1, l2, r2 =  node_indices(tree, index)
            l1_2, r1_2, l2_2, r2_2 =  node_indices(tree2, index2)
            for i in l1:r1
                idx = tree.data_reordered ? i : tree.indices[i]
                for j in l1_2:r1_2
                    idx2 = tree2.data_reordered ? j : tree2.indices[j]
                    dist_d = euclidean_distance_red(tree, idx, tree2, idx2)
                    if dist_d < r
                        push!(idx_in_ball[idx], idx2)
                    end
                end
            end
        else # tree: leaf, tree2: internal
            _in_ball(tree, index, tree2, getleft(index2), r, idx_in_ball, min_dist, max_dist)
            _in_ball(tree, index, tree2, getright(index2), r, idx_in_ball, min_dist, max_dist)
        end
    else # tree: internal
        if isleaf(tree2, index2) # tree2: leaf
            _in_ball(tree, getleft(index), tree2, index2, r, idx_in_ball, min_dist, max_dist)
            _in_ball(tree, getright(index), tree2, index2, r, idx_in_ball, min_dist, max_dist)
        else # both internal
            _in_ball(tree, getleft(index), tree2, getleft(index2), r, idx_in_ball, min_dist, max_dist)
            _in_ball(tree, getright(index), tree2, getright(index2), r, idx_in_ball, min_dist, max_dist)
            _in_ball(tree, getleft(index), tree2, getright(index2), r, idx_in_ball, min_dist, max_dist)
            _in_ball(tree, getright(index), tree2, getleft(index2), r, idx_in_ball, min_dist, max_dist)
        end
    end
    return
end

# Adds everything in this subtree since we have determined
# that the hyper rectangle completely encloses the hyper sphere
function addall(tree::KDTree, index::Int, tree2::KDTree, index2::Int, idx_in_ball::Vector{Vector{Int}})
    l1, r1, l2, r2 = node_indices(tree, index)
    l1_2, r1_2, l2_2, r2_2 = node_indices(tree2, index2)
    for r in (l1:r1, l2:r2)
        for i in r
            idx = tree.data_reordered ? i : tree.indices[i]
            for r2 in (l1_2:r1_2, l2_2:r2_2)
                for j in r2
                    idx2 = tree2.data_reordered ? j : tree2.indices[j]
                    push!(idx_in_ball[idx], idx2)
                end
            end
        end
    end
    return
end

####################################################################
# Taken from https://github.com/JuliaLang/julia/blob/v0.3.5/base/sort.jl
# and modified because I couldn't figure out how to get rid of
# the overhead when I passed in a new anonymous function
# to the "by" argument in each node. I also removed the return value.
function _select!{T <: AbstractFloat}(v::AbstractVector, k::Int, lo::Int,
                                          hi::Int, data::Matrix{T}, dim::Int)
    lo <= k <= hi || error("select index $k is out of range $lo:$hi")
     @inbounds while lo < hi
        if hi-lo == 1
            if data[dim, v[hi]] < data[dim, v[lo]]
                v[lo], v[hi] = v[hi], v[lo]
            end
            return
        end
        pivot = v[(lo+hi)>>>1]
        i, j = lo, hi
        while true
            while data[dim, v[i]] < data[dim, pivot]; i += 1; end
            while  data[dim, pivot] <  data[dim, v[j]] ; j -= 1; end
            i <= j || break
            v[i], v[j] = v[j], v[i]
            i += 1; j -= 1
        end
        if k <= j
            hi = j
        elseif i <= k
            lo = i
        else
            return
        end
    end
    return
end

# In place heap sort
function heap_sort_inplace!{T <: AbstractFloat}(xs::AbstractArray{T}, xis::AbstractArray{Int})
    @inbounds for i in length(xs):-1:2
        xs[i], xs[1] = xs[1], xs[i]
        xis[i], xis[1] = xis[1], xis[i]
        percolate_down!(xs, xis, xs[1], xis[1], i-1)
    end
end

# Binary min-heap percolate down.
function percolate_down!{T <: AbstractFloat}(xs::AbstractArray{T},
                                             xis::AbstractArray{Int},
                                             dist::T,
                                             index::Int,
                                             len::Int=length(xs))
    i = 1
    @inbounds while (l = getleft(i)) <= len
        r = getright(i)
        j = ifelse(r > len || (xs[l] > xs[r]), l, r)
        if xs[j] > dist
            xs[i] = xs[j]
            xis[i] = xis[j]
            i = j
        else
            break
        end
    end
    xs[i] = dist
    xis[i] = index
end


text:
using Symbolics, SymbolicUtils

function solve_newton(f, x, x₀; abstol=1e-8, maxiter=50)
    xₙ = Float64(x₀)
    fₙ₊₁ = x - f / Symbolics.derivative(f, x)

    for i = 1:maxiter
        xₙ₊₁ = substitute(fₙ₊₁, Dict(x => xₙ))
        if abs(xₙ₊₁ - xₙ) < abstol
            return xₙ₊₁
        else
            xₙ = xₙ₊₁
        end
    end
    return xₙ₊₁
end

n = 2
@variables ϵ a[1:n]

x = 1 + a[1]*ϵ + a[2]*ϵ^2

  eq = x^5 + ϵ*x - 1

expand(eq)

function collect_powers(eq, x, ns; max_power=100)
    eq = substitute(expand(eq), Dict(x^j => 0 for j=last(ns)+1:max_power))

    eqs = []
    for i in ns
        powers = Dict(x^j => (i==j ? 1 : 0) for j=1:last(ns))
        push!(eqs, substitute(eq, powers))
    end
    eqs
end

eqs = collect_powers(eq, ϵ, 1:2)

substitute(expand(eq), Dict(
  ϵ => 0,
  ϵ^2 => 1,
  ϵ^3 => 0,
  ϵ^4 => 0,
  ϵ^5 => 0,
  ϵ^6 => 0,
  ϵ^7 => 0,
  ϵ^8 => 0)
)

function solve_coef(eqs, ps)
    vals = Dict()

    for i = 1:length(ps)
        eq = substitute(eqs[i], vals)
        vals[ps[i]] = Symbolics.solve_for(eq ~ 0, ps[i])
    end
    vals
end

solve_coef(eqs, a)

X = 𝜀 -> 1 + a[1]*𝜀 + a[2]*𝜀^2

def_taylor(x, ps) = sum([a*x^i for (i,a) in enumerate(ps)])
def_taylor(x, ps, p₀) = p₀ + def_taylor(x, ps)

n = 3
@variables ϵ M a[1:n]
x = def_taylor(ϵ, a, M)

expand_sin(x, n) = sum([(isodd(k) ? -1 : 1)*(-x)^(2k-1)/factorial(2k-1) for k=1:n])

expand_sin(0.1, 10) ≈ sin(0.1)

eq = x - ϵ * expand_sin(x, n) - M

eqs = collect_powers(eq, ϵ, 1:n)

vals = solve_coef(eqs, a)

x′ = substitute(x, vals)
X = (𝜀, 𝑀) -> substitute(x′, Dict(ϵ => 𝜀, M => 𝑀))
X(0.01671, π/2)

ext: jl


text:
"""
    pairwise(dist::StringDistance, xs::AbstractVector, ys::AbstractVector = xs; preprocess = true)

Compute distances between all pairs of elements in `xs`  and `ys` according to the
`StringDistance` `dist`. Returns a matrix R such that `R[i, j]` corrresponds to the distance between `xs[i]` and `ys[j]`.

Set `preprocess` to false if no preprocessing should be used.

Both symmetric and asymmetric versions are available.

### Examples
```julia-repl
julia> using StringDistances
julia> iter = ["New York", "Princeton"]
julia> pairwise(Levenshtein(), iter)
2×2 Array{Float64,2}:
 0.0  9.0
 9.0  0.0
julia> iter2 = ["San Francisco"]
julia> pairwise(Levenshtein(), iter, iter2)
2×1 Array{Float64,2}:
 12.0
 10.0
```
"""
function StatsAPI.pairwise(dist::Union{StringSemiMetric, StringMetric}, xs::AbstractVector, ys::AbstractVector = xs; preprocess = true)
    T = result_type(dist, eltype(xs), eltype(ys))
    R = Matrix{T}(undef, length(xs), length(ys))
    pairwise!(R, dist, xs, ys; preprocess = preprocess)
end

"""
    pairwise!(R::AbstractMatrix, dist::Union{StringSemiMetric, StringMetric}, xs::AbstractVector, ys::AbstractVector = xs; preprocess = true)

Compute distances between all pairs of elements in `xs` and `ys` according to the
`Union{StringSemiMetric, StringMetric}` `dist` and write the result in `R`. `R[i, j]` corresponds to the distance between `xs[i]` and `ys[j]`.

Set `preprocess` to false if no preprocessing should be used.
"""
function StatsAPI.pairwise!(R::AbstractMatrix, dist::Union{StringSemiMetric, StringMetric}, xs::AbstractVector, ys::AbstractVector = xs; preprocess = true)
    length(xs) == size(R, 1) || throw(DimensionMismatch("inconsistent length"))
    length(ys) == size(R, 2) || throw(DimensionMismatch("inconsistent length"))
    (xs === ys) ?
        _symmetric_pairwise!(R, dist, xs; preprocess = preprocess) :
        _asymmetric_pairwise!(R, dist, xs, ys; preprocess = preprocess)
end

function _symmetric_pairwise!(R::AbstractMatrix, dist::Union{StringSemiMetric, StringMetric}, xs::AbstractVector; preprocess = true)
    if preprocess
        xs = _preprocess_list(dist, xs)
    end
    for i in 1:length(xs)
        # handle missing
        R[i, i] = xs[i] != xs[i]
        Threads.@threads for j in (i+1):length(xs)
            R[i, j] = R[j, i] = evaluate(dist, xs[i], xs[j])
        end
    end
    return R
end

function _asymmetric_pairwise!(R::AbstractMatrix, dist::Union{StringSemiMetric, StringMetric}, xs::AbstractVector, ys::AbstractVector; preprocess = true)
    if preprocess
        objxs = _preprocess_list(dist, xs)
        objys = xs === ys ? objxs : _preprocess_list(dist, ys)
    else
        objxs = xs
        objys = ys
    end
    for i in 1:length(objxs)
        Threads.@threads for j in 1:length(objys)
            R[i, j] = evaluate(dist, objxs[i], objys[j])
        end
    end
    return R
end

_preprocess_list(dist::Union{StringSemiMetric, StringMetric}, xs)  = xs
_preprocess_list(dist::AbstractQGramDistance, xs) = fetch.(map(x -> (Threads.@spawn x === missing ? x : QGramSortedVector(x, dist.q)), xs))

max_stars_repo_name: matthieugomez/StringDistances

max_issues_repo_name: matthieugomez/StringDistances

max_forks_repo_name: matthieugomez/StringDistances


text:
#
# Linear Algebra of order 2 Tensors
#
# Even though DenseTensor{_,2} is strided
# and passable to BLAS/LAPACK, it cannot
# be made <: StridedArray

function Base.:*(T1::Tensor{ElT1,2,StoreT1,IndsT1},
                 T2::Tensor{ElT2,2,StoreT2,IndsT2}) where
                                       {ElT1,StoreT1<:Dense,IndsT1,
                                        ElT2,StoreT2<:Dense,IndsT2}
  RM = matrix(T1)*matrix(T2)
  indsR = IndsT1(ind(T1,1),ind(T2,2))
  return Tensor(Dense{promote_type(ElT1,ElT2)}(vec(RM)),indsR)
end

function LinearAlgebra.exp(T::DenseTensor{ElT,2}) where {ElT}
  expTM = exp(matrix(T))
  return Tensor(Dense{ElT}(vec(expTM)),inds(T))
end

function expHermitian(T::DenseTensor{ElT,2}) where {ElT}
  # exp(::Hermitian/Symmetric) returns Hermitian/Symmetric,
  # so extract the parent matrix
  expTM = parent(exp(Hermitian(matrix(T))))
  return Tensor(Dense{ElT}(vec(expTM)),inds(T))
end

# svd of an order-2 tensor
function LinearAlgebra.svd(T::DenseTensor{ElT,2,IndsT};
                           kwargs...) where {ElT,IndsT}
  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  fastSVD::Bool = get(kwargs,:fastSVD,false)

  if fastSVD
    MU,MS,MV = svd(matrix(T))
  else
    MU,MS,MV = recursiveSVD(matrix(T))
  end
  conj!(MV)

  P = MS.^2
  truncate!(P;mindim=mindim,
              maxdim=maxdim,
              cutoff=cutoff,
              absoluteCutoff=absoluteCutoff,
              doRelCutoff=doRelCutoff)
  dS = length(P)
  if dS < length(MS)
    MU = MU[:,1:dS]
    resize!(MS,dS)
    MV = MV[:,1:dS]
  end

  # Make the new indices to go onto U and V
  u = eltype(IndsT)(dS)
  v = eltype(IndsT)(dS)
  Uinds = IndsT((ind(T,1),u))
  Sinds = IndsT((u,v))
  Vinds = IndsT((ind(T,2),v))
  U = Tensor(Dense{ElT}(vec(MU)),Uinds)
  S = Tensor(Diag{Vector{real(ElT)}}(MS),Sinds)
  V = Tensor(Dense{ElT}(vec(MV)),Vinds)
  return U,S,V
end

function eigenHermitian(T::DenseTensor{ElT,2,IndsT};
                        kwargs...) where {ElT,IndsT}
  ispossemidef::Bool = get(kwargs,:ispossemidef,false)
  maxdim::Int = get(kwargs,:maxdim,minimum(dims(T)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)

  DM,UM = eigen(Hermitian(matrix(T)))

  # Sort by largest to smallest eigenvalues
  p = sortperm(DM; rev = true)
  DM = DM[p]
  UM = UM[:,p]

  if ispossemidef
    truncate!(DM;maxdim=maxdim,
                 cutoff=cutoff,
                 absoluteCutoff=absoluteCutoff,
                 doRelCutoff=doRelCutoff)
    dD = length(DM)
    if dD < size(UM,2)
      UM = UM[:,1:dD]
    end
  else
    dD = length(DM)
  end

  # Make the new indices to go onto U and V
  u = eltype(IndsT)(dD)
  v = eltype(IndsT)(dD)
  Uinds = IndsT((ind(T,1),u))
  Dinds = IndsT((u,v))
  U = Tensor(Dense{ElT}(vec(UM)),Uinds)
  D = Tensor(Diag{Vector{real(ElT)}}(DM),Dinds)
  return U,D
end

function LinearAlgebra.qr(T::DenseTensor{ElT,2,IndsT}) where {ElT,
                                                              IndsT}
  # TODO: just call qr on T directly (make sure
  # that is fast)
  QM,RM = qr(matrix(T))
  # Make the new indices to go onto Q and R
  q,r = inds(T)
  q = dim(q) < dim(r) ? sim(q) : sim(r)
  Qinds = IndsT((ind(T,1),q))
  Rinds = IndsT((q,ind(T,2)))
  Q = Tensor(Dense{ElT}(vec(Matrix(QM))),Qinds)
  R = Tensor(Dense{ElT}(vec(RM)),Rinds)
  return Q,R
end

function polar(T::DenseTensor{ElT,2,IndsT}) where {ElT,IndsT}
  QM,RM = polar(matrix(T))
  dim = size(QM,2)
  # Make the new indices to go onto Q and R
  q = eltype(IndsT)(dim)
  # TODO: use push/pushfirst instead of a constructor
  # call here
  Qinds = IndsT((ind(T,1),q))
  Rinds = IndsT((q,ind(T,2)))
  Q = Tensor(Dense{ElT}(vec(QM)),Qinds)
  R = Tensor(Dense{ElT}(vec(RM)),Rinds)
  return Q,R
end


text:
struct Solution
    ret_val::Int
    primal::Vector{Float64}
    dual_eq::Vector{Float64}
    dual_ineq::Vector{Float64}
    slack::Vector{Float64}
    objval::Float64
end
Solution() = Solution(0,
                      Float64[], Float64[], Float64[], Float64[], NaN)

mutable struct Data
    m::Int
    n::Int
    IA::Vector{Int}
    JA::Vector{Int}
    VA::Vector{Float64}
    b::Vector{Float64}
    IG::Vector{Int}
    JG::Vector{Int}
    VG::Vector{Float64}
    h::Vector{Float64}
    objconstant::Float64
    c::Vector{Float64}
end

mutable struct Cone
    f::Int # number of linear equality constraints
    l::Int # length of LP cone
    q::Int # length of SOC cone
    qa::Vector{Int} # array of second-order cone constraints
    ep::Int # number of primal exponential cone triples
    eqsetconstant::Dict{Int, Float64} # For the constant of EqualTo
    eqnrows::Dict{Int, Int} # The number of rows of Zeros
    ineqsetconstant::Dict{Int, Float64} # For the constant of LessThan and GreaterThan
    ineqnrows::Dict{Int, Int} # The number of rows of each vector sets except Zeros
    function Cone()
        new(0, 0, 0, Int[], 0,
            Dict{Int, Float64}(),
            Dict{Int, UnitRange{Int}}(),
            Dict{Int, Float64}(),
            Dict{Int, UnitRange{Int}}())
    end
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.3'
#       jupytext_version: 0.8.6
#   kernelspec:
#     display_name: Julia 1.0.3
#     language: julia
#     name: julia-1.0
# ---

module DiscreteTimeSEIRD
using DifferentialEquations

function seird(du, u, parms, t)
    S,E,Ir,Id,R,D = u
    β,δ,γ,Γ,μ,ϵ,ω = parms
    I = Ir + Id
    N = S + E + I + R + D
    dS = S - β*S*I/N + ω*R
    dE = E + β*S*I/N - γ*E + ϵ
    dIr = Ir + γ*(1-μ)E - γ*Ir + ϵ
    dId = Id + γ*μ*E - Γ*Id + ϵ
    dR = R + γ*Ir - ω*R
    dD = D + Γ*Id
    du = [dS, dE, dIr, dId, dR, dD]
end

end


text:
using NNlib, Test

@testset "NNlib" begin

include("activation.jl")
include("conv.jl")

xs = [-100_000, -100_000.]
@test softmax(xs) ≈ [0.5, 0.5]
@test logsoftmax(xs) ≈ log.([0.5, 0.5])

xs = rand(5)
@test softmax(xs) ≈ exp.(xs) ./ sum(exp.(xs))
@test logsoftmax(xs) ≈ log.(softmax(xs))
@test logsigmoid.(xs) ≈ log.(sigmoid.(xs))

xs = rand(5,10)
@test softmax(xs) ≈ exp.(xs) ./ sum(exp.(xs), dims = 1)
@test logsoftmax(xs) ≈ log.(softmax(xs))
@test logsigmoid.(xs) ≈ log.(sigmoid.(xs))

for T in [:Float32, :Float64]
  @eval @test logsigmoid.($T[-100_000, 100_000.]) ≈ $T[-100_000, 0.]
end

## compare the outputs with the PyTorch nn.LogSoftmax returns
xs = Float32[1, 2, 3000.]
@test logsoftmax(xs) ≈ [-2999, -2998, 0]

xs = Float32[1 2 3; 1000 2000 3000]
@test logsoftmax(xs) ≈ [-999 -1998 -2997; 0 0 0.]

@test NNlib.∇logsoftmax(ones(size(xs)), xs) ≈ zeros(Float32, size(xs))
@test NNlib.∇softmax(ones(size(xs)), xs) ≈ zeros(Float32, size(xs))

xs = [-0.238639 0.748142 -0.283194 -0.525461 -1.5348 -0.797842; 0.690384 0.211427 0.254794 -0.213572 -0.314174 -0.372663; -1.14637 -0.577988 0.718952 0.91972 -0.620773 0.929977]
@test isapprox(NNlib.∇logsoftmax(ones(size(xs)), xs), [0.237703 -0.621474 0.448193 0.546047 0.564185 0.632273; -0.930163 0.0519798 0.0549979 0.3799 -0.477112 0.437428; 0.69246 0.569494 -0.503191 -0.925947 -0.0870738 -1.0697]; rtol = 1e-6)
@test isapprox(NNlib.∇softmax(ones(size(xs)), xs), zeros(size(xs)); atol = 1e-6)
end


text:
module Transducers

export ttake, tmap, tfilter, partition_all, random_sample, transduce,
       dedupe, treplace, tpush!

type Reduced
  val
end

function ensure_reduced(x::Reduced)
  x
end

function ensure_reduced(x)
  Reduced(x)
end

function unreduced(x::Reduced)
  x.val
end

function unreduced(x)
  x
end

function treduce(fn, iterable, initializer=None)
  if initializer == None
    accum_value = fn()
  else
    accum_value = initializer
  end

  for x in iterable
    accum_value = fn(accum_value, x)
    if typeof(accum_value) == Reduced
      return accum_value.val
    end
  end
  accum_value
end

function tpush!()
  Any[]
end

function tpush!(r)
  r
end

function tpush!(r, x)
  push!(r, x)
end

"""
Transduce entry point - reduce collection into start using
reducing function f, applying xform with each step.
"""
function transduce(xform, f, start)
  transduce(xform, f, f(), start)
end

function transduce(xform, f, start, coll)
  let reducer = xform(f)
    reducer((treduce(reducer, coll, start)))
  end
end

function tmap(f)
  function _map_xducer(step)
    _map_step() = step()
    _map_step(r) = step(r)
    _map_step(r, x) = step(r, f(x))
    _map_step
  end
  _map_xducer
end

function tfilter(pred)
  function _filter_xducer(step)
    _filter_step() = step()
    _filter_step(r) = step(r)
    _filter_step(r, x) = pred(x) ? step(r, x) : r
    _filter_step
  end
  _filter_xducer
end

function ttake(n::Int)
  function _take_xducer(step)
    counter = n
    _take_step() = step()
    _take_step(r) = step(r)
    function _take_step(r, x)
      old_n = counter
      counter -= 1
      r = old_n > 0 ? step(r, x) : r # ensure_reduced
      counter <= 0 ? ensure_reduced(r) : r
    end
    _take_step
  end
  _take_xducer
end

function partition_all(n::Int)
  function _partition_all_xducer(step)
    temp = Any[]
    _partition_all_step() = step()

    function _partition_all_step(r)
      if size(temp)[1] == 0
        r
      else
        _temp = temp[:] #benchmark vs. copy?
        empty!(temp)
        _r =  unreduced(step(r, _temp))
        step(_r)
      end
    end

    function _partition_all_step(r, x)
      push!(temp, x)
      if size(temp)[1] == n
        _temp = temp[:]
        empty!(temp)
        step(r, _temp)
      else
        r
      end
    end

    _partition_all_step
  end
  _partition_all_xducer
end

function treplace(smap::Dict{Any,Any})
  function _replace_xducer(step)
    _replace_step() = step()
    _replace_step(r) = step(r)
    _replace_step(r, x) = haskey(smap, x) ? step(r, smap[x]) : step(r, x)
    _replace_step
  end
  _replace_xducer
end

function dedupe(step)
  prev = None
  _dedupe_step() = step()
  _dedupe_step(r) = step(r)
  function _dedupe_step(r, x)
    if prev == None
      prev = x
      step(r, x)
    elseif x != prev
      prev = x
      step(r, x)
    else
      r
    end
  end
  _dedupe_step
end

function random_sample(prob::FloatingPoint)
  function _random_sample_xducer(step)
    _random_sample_step() = step()
    _random_sample_step(r) = step(r)
    _random_sample_step(r, x) = rand() < prob ? step(r, x)  : r
    _random_sample_step
  end
  _random_sample_xducer
end

# module end
end


text:
# show

function show_readable(io::IO, x::T) where {T <: Real}
    str = readable(x)
    print(io, str)
end

function show_readable(x::T) where {T <: Real}
    str = readable(x)
    print(stdout, str)
end


text:
using Oceananigans.Utils: with_tracers
using Oceananigans.TurbulenceClosures: with_tracers

"""
    ModelForcing(; u=zeroforcing, v=zeroforcing, w=zeroforcing, tracer_forcings...)

Return a named tuple of forcing functions for each solution field.

Example
=======

julia> u_forcing = SimpleForcing((x, y, z, t) -> exp(z) * cos(t))

julia> model = Model(forcing=ModelForcing(u=u_forcing))
"""
function ModelForcing(; u=zeroforcing, v=zeroforcing, w=zeroforcing, tracer_forcings...)
    u = at_location((Face, Cell, Cell), u)
    v = at_location((Cell, Face, Cell), v)
    w = at_location((Cell, Cell, Face), w)

    return merge((u=u, v=v, w=w), tracer_forcings)
end

at_location(location, u::Function) = u
at_location(location, u::SimpleForcing) =
    SimpleForcing{location[1], location[2], location[3]}(u.func, u.parameters)

default_tracer_forcing(args...) = zeroforcing
ModelForcing(tracers, proposal_forcing) = with_tracers(tracers, proposal_forcing, default_tracer_forcing,
                                                       with_velocities=true)


text:
function factorial4(n)
  if n == 0
    return BigInt(1)
  else
    return n * factorial4(n-1)
  end
end

nCr :: Int16 -> Int16 -> BigInt
function nCr(n, r)
  return BigInt(factorial4(n)/(factorial4(n-r) * factorial4(r)))
end

get_draw :: Int16 -> Float64 -> Int16 -> Float64
get_draw(n, p, x) = Float64(nCr(n, x) * p^x * (1-p)^(n-x))
println(nCr(100, 10))
println(get_draw(5, 1/2, 1))

function binomial_rv2(n, p)
  m = 0
  for i in range(1, n)
    if p < rand()
      m += 1
    end
  end
  return m
end

function binomial_rv3(n, p)
  rands = [rand() for i in range(1, n)]
  return length(rands[rands.<p])
end

binomial_rv4 = (n, p) -> length(filter(x -> x < p, [rand() for i in range(1, n)]))
#binomial_rv3 = (n, p) -> sum(filter(_ -> 1/2 < rand(), [1]*n))
#println(binomial(5, 2))
#println(filter(prob(1/2), [1, 1]))
#println(binomial_rv2(5, 1/2))
#println(binomial_rv3(5, 1/2))
println(binomial_rv4(5, 1/2))


text:
module _Taylor

include("euler_taylor_simple.jl")
include("euler_taylor_inplace.jl")

end


text:
function Base.tryparse(::Type{InfExtendedReal{T}}, str::AbstractString) where T
  val = tryparse(Infinite, str)
  if val === nothing
    val = tryparse(T, str)
  end
  return val !== nothing ? InfExtendedReal{T}(val) : nothing
end


text:
type DeckDesign{T<:Integer,U<:String}
    rlen::T
    slen::T
    ranks::Array{U,1}
    suits::Array{U,1}
    hlen::T
end

type Deck{T<:Integer}
    cards::Array{T,1}
    design::DeckDesign
end

Deck(n::Integer, des::DeckDesign) = Deck([n], des)

function pokerlayout()
    r = [map(string, 2:10), "J", "Q", "K", "A"]
    r = map(utf8, r)
    s = ["\u2663", "\u2666", "\u2665", "\u2660"]
    DeckDesign(13, 4, r, s, 5)
end

function fresh(des::DeckDesign)
    Deck(collect(1:des.rlen*des.slen), des)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
##############################
#
# obj-Files
#
##############################

function load(io::Stream{format"OBJ"}, MeshType::Type{MT} = GLNormalMesh) where {MT <: AbstractMesh}
    Tv,Tn,Tuv,Tf = vertextype(MT), normaltype(MT), texturecoordinatetype(MT), facetype(MT)
    v,n,uv,f     = Tv[], Tn[], Tuv[], Tf[]
    f_uv_n_faces = (f, GLTriangle[], GLTriangle[])
    last_command = ""
    attrib_type  = nothing
    for full_line::String in eachline(stream(io))
        # read a line, remove newline and leading/trailing whitespaces
        line = strip(chomp(full_line))
        !isascii(line) && error("non valid ascii in obj")

        if !startswith(line, "#") && !isempty(line) && !all(iscntrl, line) #ignore comments
            lines   = split(line)
            command = popfirst!(lines) #first is the command, rest the data

            if "v" == command # mesh always has vertices
                push!(v, Point{3, Float32}(parse.(Float32, lines))) # should automatically convert to the right type in vertices(mesh)
            elseif "vn" == command && hasnormals(MT)
                push!(n, Normal{3, Float32}(parse.(Float32, lines)))
            elseif "vt" == command && hastexturecoordinates(MT)
                if length(lines) == 2
                    push!(uv, UV{Float32}(parse.(Float32, lines)))
                elseif length(lines) == 3
                    push!(uv, UVW{Float32}(parse.(Float32, lines)))
                else
                    error("Unknown UVW coordinate: $lines")
                end
            elseif "f" == command #mesh always has faces
                if any(x->occursin("//", x), lines)
                    fs = process_face_normal(lines)
                elseif any(x->occursin("/", x), lines)
                    fs = process_face_uv_or_normal(lines)
                else
                    append!(f, triangulated_faces(Tf, lines))
                    continue
                end
                for i = 1:length(first(fs))
                    append!(f_uv_n_faces[i], triangulated_faces(Tf, getindex.(fs, i)))
                end
            else
                #TODO
            end
        end
    end
    attributes = Dict{Symbol, Any}()
    !isempty(f) && (attributes[:faces] = f)
    !isempty(v) && (attributes[:vertices] = v)
    if !isempty(n)
        attributes[:normals] = if !isempty(f_uv_n_faces[3])
            _n = similar(v, eltype(n))
            for (vf, nf) in zip(f, f_uv_n_faces[3])
                for (vidx, nidx) in zip(vf, nf)
                    _n[vidx] = n[nidx]
                end
            end
            _n
        else
            # these are not per vertex normals, which we
            # can't deal with at the moment
            if length(v) != length(n)
                normals(v, f, Tn)
            else
                n
            end
        end
    end
    if !isempty(uv)
        attributes[:texturecoordinates] = if !isempty(f_uv_n_faces[2])
            _uv = similar(v, eltype(uv))
            for (vf, uvf) in zip(f, f_uv_n_faces[2])
                for (vidx, uvidx) in zip(vf, uvf)
                    _uv[vidx] = uv[uvidx]
                end
            end
            _uv
        else
            uv
        end
    end
    return MT(GeometryTypes.homogenousmesh(attributes))::MT
end

# of form "f v1 v2 v3 ....""
process_face(lines::Vector{S}) where {S <: AbstractString} = (lines,) # just put it in the same format as the others
# of form "f v1//vn1 v2//vn2 v3//vn3 ..."
process_face_normal(lines::Vector{S}) where {S <: AbstractString} = split.(lines, "//")
# of form "f v1/vt1 v2/vt2 v3/vt3 ..." or of form "f v1/vt1/vn1 v2/vt2/vn2 v3/vt3/vn3 ...."
process_face_uv_or_normal(lines::Vector{S}) where {S <: AbstractString} = split.(lines, Ref('/'))

function triangulated_faces(::Type{Tf}, vertex_indices::Vector{<:AbstractString}) where {Tf}
    poly_face = Face{length(vertex_indices), UInt32}(parse.(UInt32, vertex_indices))
    decompose(Tf, poly_face)
end


text:
"""
    conductance(h::Hypergraph, subset::Set{Int})::Float64

Calculate unweighted hypergraph conductance of `subset`.
note: ∅ ⊊ `subset` ⊊ `1:nhv(h)`

For more information see `1. Introduction` at:
Spectral Properties of Hypergraph Laplacian and Approximation Algorithms
auhtors: T-H. Hubert Chan, Anand Louis, Zhihao Gavin Tang, and Chenzi Zhang
"""
function conductance(h::Hypergraph, subset::Set{Int})::Float64
  if isempty(subset) error("`subset` is not allowed empty.") end
  if subset == Set(1:nhv(h)) error("`subset` is not allowed true subset of `h`.") end
  subset2 = setdiff(Set(1:nhv(h)), subset)
  volS = sum([sum(values(gethyperedges(h, node))) for node in subset])
  volV_S = sum([sum(values(gethyperedges(h, node))) for node in subset2])
  cutS = 0
  for he in 1:nhe(h)
    he_vertices = keys(getvertices(h, he))
    if isempty(intersect(he_vertices, subset)) continue end
    if isempty(intersect(he_vertices, subset2)) continue end
    cutS += 1
  end
  return cutS / min(volS, volV_S)
end


text:
struct gpuTrackingResults{
        TS <: Union{gpuTrackingState, sgpuTrackingState},
        C <: AbstractCorrelator
    }
    state::TS
    correlator::C
    got_correlator::Bool
    bit_buffer::BitBuffer
    cn0::typeof(1.0dBHz)
end

@inline get_state(results::gpuTrackingResults) = results.state

@inline get_carrier_doppler(results::gpuTrackingResults) = get_carrier_doppler(results.state)

@inline get_carrier_phase(results::gpuTrackingResults) = get_carrier_phase(results.state) * 2π

@inline get_code_doppler(results::gpuTrackingResults) = get_code_doppler(results.state)

@inline get_code_phase(results::gpuTrackingResults) = get_code_phase(results.state)

@inline get_correlator(results::gpuTrackingResults) = results.correlator

@inline get_early(results::gpuTrackingResults) = get_early(get_correlator(results))

@inline get_prompt(results::gpuTrackingResults) = get_prompt(get_correlator(results))

@inline get_late(results::gpuTrackingResults) = get_late(get_correlator(results))

@inline get_bits(results::gpuTrackingResults) = get_bits(results.bit_buffer)

@inline get_num_bits(results::gpuTrackingResults) = length(results.bit_buffer)

@inline get_cn0(results::gpuTrackingResults) = results.cn0

@inline get_secondary_code_or_bit_found(results::TrackingResults) =
    found(get_sc_bit_detector(results.state))


text:
#
# ITensors.jl extensions
#

# Generalize siteind to n-dimensional lattice
function ITensors.siteind(st::SiteType, N1::Integer, N2::Integer, Ns::Integer...; kwargs...)
  s = siteind(st; kwargs...)
  if !isnothing(s)
    ts = "n1=$N1,n2=$N2"
    for i in eachindex(Ns)
      ts *= ",n$(i + 2)=$(Ns[i])"
    end
    return addtags(s, ts)
  end
  isnothing(s) && error(space_error_message(st))
end

# Generalize siteinds to n-dimensional lattice
function ITensors.siteinds(str::AbstractString, N1::Integer, N2::Integer, Ns::Integer...; kwargs...)
  st = SiteType(str)
  return [siteind(st, ns...) for ns in Base.product(1:N1, 1:N2, UnitRange.(1, Ns)...)]
end

#
# HyperCubic lattice
#

# An N-dimensional hypercubic lattice with periodic
# boundary conditions
struct HyperCubic{N}
  dims::NTuple{N,Int}
end
Base.size(l::HyperCubic) = l.dims

const Chain = HyperCubic{1}
const Square = HyperCubic{2}
const Cubic = HyperCubic{3}

struct Edge{N}
  edge::Tuple{NTuple{N,Int},NTuple{N,Int}}
  boundary::Bool
end

Base.getindex(edge::Edge, n::Int) = edge.edge[n]

Base.reverse(edge::Edge) = Edge(reverse(edge.edge), edge.boundary)

sites(lattice::HyperCubic) = (Tuple(s) for s in CartesianIndices(axes(lattice)))

function coordinate_tag(n)
  str = replace("$n", ")" => "")
  str = replace(str, "(" => "")
  str = replace(str, " " => "")
  if length(n) > 1
    str = replace(str, "," => ".")
  else
    str = replace(str, "," => "")
  end
  return str
end

function link_tag(n1, n2)
  link_string = "$(coordinate_tag(n1))↔$(coordinate_tag(n2))"
  start_ind = nextind(link_string, 0, 1)
  stop_ind = min(ncodeunits(link_string), nextind(link_string, 0, 8))
  link_string = link_string[start_ind:stop_ind]
  return TagSet(link_string)
end

function ITensors.linkinds(lattice::HyperCubic; linkdims, addtags=ts"")
  dims = size(lattice)
  N = length(dims)
  linkinds_dict = Dict{Edge{N},Index{typeof(linkdims)}}()
  for n in sites(lattice), edge_n in incident_edges(lattice, n)
    l = Index(linkdims; tags=ITensors.addtags(link_tag(edge_n.edge...), addtags))
    get!(linkinds_dict, edge_n, l)
  end
  return linkinds_dict
end

using ITensors: data

function Base.isless(i1::Index, i2::Index)
  return isless((id(i1), plev(i1), tags(i1), dir(i1)), (id(i2), plev(i2), tags(i2), dir(i2)))
end

Base.isless(ts1::TagSet, ts2::TagSet) = isless(data(ts1), data(ts2))

function is_in_edge(site::Tuple, edge::Edge)
  if site == edge[1]
    return false
  elseif site == edge[2]
    return true
  else
    error("Site $site is not incident to edge $edge")
  end
end

function get_link_ind(linkinds_dict::Dict, edge::Edge, site::Tuple)
  l = linkinds_dict[edge]
  return is_in_edge(site, edge) ? dag(l) : l
end

# A network of link indices for a HyperCubic lattice, with
# no site indices.
function inds_network(dims::Int...; linkdims, kwargs...)
  site_inds = fill(Index{typeof(linkdims)}[], dims)
  return inds_network(site_inds; linkdims=linkdims, kwargs...)
end

function inds_network(site_inds::Array{<:Index,N}; kwargs...) where {N}
  return inds_network(map(x -> [x], site_inds); kwargs...)
end

# A network of link indices for a HyperCubic lattice, with
# site indices specified.
function inds_network(site_inds::Array{<:Vector{<:Index},N}; linkdims, addtags=ts"") where {N}
  dims = size(site_inds)
  lattice = HyperCubic(dims)
  linkinds_dict = linkinds(lattice; linkdims, addtags)
  inds = Array{Vector{Index{typeof(linkdims)}},N}(undef, dims)
  for n in sites(lattice)
    inds_n = [get_link_ind(linkinds_dict, edge_n, n) for edge_n in incident_edges(lattice, n)]
    inds[n...] = append!(inds_n, site_inds[n...])
  end
  return inds
end

function itensor_network(dims::Int...; linkdims)
  return ITensor.(inds_network(dims...; linkdims))
end

function onehot_tuple(n::Integer, length::Val{N}) where {N}
  return ntuple(i -> i == n ? 1 : 0, Val(N))
end

# Obtain the neighbor in dimension `dim` in direction `dir`, for example:
# neighbor((2, 2), 2, 1) == (2, 3)
# neighbor((2, 2), 2, -1) == (2, 1)
# neighbor((2, 2), 1, -1) == (1, 2)
# neighbor((2, 3), 2, 1, lattice_size=(3, 3)) == (2, 1)
function neighbor(site::NTuple{N,Int}, dim::Int, dir::Int; lattice_size=typemax(eltype(site))) where {N}
  return map((a, b, d) -> mod1(a + dir * b, d), site, onehot_tuple(dim, Val(N)), lattice_size)
end

#function isboundary(site::dir::Int,
# Check if the edge connecting to the specified neighbor of the site
# crosses the boundary of the lattice
function isboundary(site::NTuple{N,Int}, dim::Int, dir::Int; lattice_size=typemax(eltype(site))) where {N}
  return ((dir == 1) && (site[dim] == lattice_size[dim])) || ((dir == -1) && site[dim] == 1)
end

# The neighboring sites of the specified site
function neighbors(lattice::HyperCubic{N}, site::NTuple{N,Int}) where {N}
  lattice_size = size(lattice)
  site_neighbors = Vector{NTuple{N,Int}}()
  for dim in 1:N, dir in (-1, 1)
    site_neighbor = neighbor(site, dim, dir; lattice_size=lattice_size)
    push!(site_neighbors, neighbor)
  end
  return site_neighbors
end

# All of the edges connected to the vertex `site`
function incident_edges(lattice::HyperCubic{N}, site::NTuple{N,Int}) where {N}
  lattice_size = size(lattice)
  site_edges = Vector{Edge{N}}()
  for dim in 1:N, dir in (-1, 1)
    site_neighbor = neighbor(site, dim, dir; lattice_size=lattice_size)
    boundary = false
    if isboundary(site, dim, dir; lattice_size=lattice_size)
      boundary = true
    end
    edge = (site, site_neighbor)
    if dir == -1
      edge = reverse(edge)
    end
    push!(site_edges, Edge(edge, boundary))
  end
  return site_edges
end

bonds(lattice::HyperCubic) = ((s, n) for s in sites(lattice) for n in directed_neighbors(lattice, s))

function unique_bonds(lattice::HyperCubic)
  # Turn into Sets so order doesn't matter
  bond_sets = [Set(b) for b in bonds(l)]
  unique!(bond_sets)
  bond_tuples = map(Tuple, bond_sets)
  sort!(bond_tuples)
  return bond_tuples
end

#
# Contraction methods
#

# Return the right and bottom boundary projectors onto
# the local state `state`
function boundary_projectors(tn::Matrix{ITensor}, state=1)
  top_row = tn[1, :]
  bottom_row = tn[end, :]
  left_column = tn[:, 1]
  right_column = tn[:, end]
  bottom_boundary_inds = commonind.(bottom_row, top_row)
  right_boundary_inds = commonind.(right_column, left_column)
  ψr = ITensors.state.(right_boundary_inds, state)
  ψb = ITensors.state.(bottom_boundary_inds, state)
  return ψr, ψb
end

function project_boundary(tn::Matrix{ITensor}, state=1)
  Nx, Ny = size(tn)
  ψr, ψb = boundary_projectors(tn, state)
  for n in 1:Nx
    tn[n, 1] = tn[n, 1] * ψr[n]
    tn[n, end] = tn[n, end] * dag(ψr[n])
  end
  for n in 1:Ny
    tn[1, n] = tn[1, n] * ψb[n]
    tn[end, n] = tn[end, n] * dag(ψb[n])
  end
  return tn
end

Base.keytype(m::MPS) = keytype(data(m))

indtype(i::Index) = typeof(i)
indtype(T::Type{<:Index}) = T
indtype(is::Tuple{Vararg{<:Index}}) = eltype(is)
indtype(is::Vector{<:Index}) = eltype(is)
indtype(A::ITensor...) = indtype(inds.(A))

indtype(tn1, tn2) = promote_type(indtype(tn1), indtype(tn2))
indtype(tn) = mapreduce(indtype, promote_type, tn)

function filter_alllinkinds(f, tn)
  linkinds = Dict{Tuple{keytype(tn),keytype(tn)},Vector{indtype(tn)}}()
  for n in keys(tn), m in keys(tn)
   if f(n, m)
     is = commoninds(tn[n], tn[m])
     if !isempty(is)
       linkinds[(n, m)] = is
     end
   end
 end
 return linkinds
end
alllinkinds(tn) = filter_alllinkinds(≠, tn)
inlinkinds(tn) = filter_alllinkinds(>, tn)
outlinkinds(tn) = filter_alllinkinds(<, tn)

function split_links(H::Union{MPS,MPO}; split_tags=("" => ""), split_plevs=(0 => 1))
  left_tags, right_tags = split_tags
  left_plev, right_plev = split_plevs
  l = outlinkinds(H)
  Hsplit = copy(H)
  for bond in keys(l)
    n1, n2 = bond
    lₙ = l[bond]
    left_l_n = setprime(addtags(lₙ, left_tags), left_plev)
    right_l_n = setprime(addtags(lₙ, right_tags), right_plev)
    Hsplit[n1] = replaceinds(Hsplit[n1], lₙ => left_l_n)
    Hsplit[n2] = replaceinds(Hsplit[n2], lₙ => right_l_n)
  end
  return Hsplit
end

split_links(H::Vector{ITensor}, args...; kwargs...) = data(split_links(MPS(H), args...; kwargs...))

function insert_projectors(H, psi, projs)
  N = length(H)
  Hpsi = Vector{ITensor}(undef, N)
  Hpsi[1] = psi[1] * H[1] * projs[1][1]
  for n in 2:(N - 1)
    Hpsi[n] = projs[n - 1][2] * psi[n] * H[n] * projs[n][1]
  end
  Hpsi[N] = projs[N - 1][2] * psi[N] * H[N]
  return Hpsi
end

# Compute the truncation projectors for the network,
# contracting from top to bottom
function truncation_projectors(tn::Matrix{ITensor}; maxdim=maxdim_arg(tn), cutoff=1e-8, split_tags=("" => ""), split_plevs=(0 => 1))
  nrows, ncols = size(tn)
  U = Matrix{ITensor}(undef, nrows - 3, ncols - 1)
  Ud = copy(U)
  tn_split = copy(tn)
  if nrows ≤ 3
    return tn_split, U, Ud
  end
  psi = tn[1, :]
  psi_split = split_links(psi; split_tags=split_tags, split_plevs=split_plevs)
  tn_split[1, :] .= psi_split
  for n in 1:(nrows - 3)
    H = tn[n + 1, :]
    projs = truncation_projectors(H, psi; split_tags=split_tags, split_plevs=split_plevs, cutoff=cutoff)
    U[n, :] = first.(projs)
    Ud[n, :] = last.(projs)
    if n < size(U, 1)
      for m in 1:size(U, 2)
        U[n, m], Ud[n, m] = split_links([U[n, m], Ud[n, m]]; split_tags=split_tags, split_plevs=split_plevs)
      end
    end
    H_split = split_links(H; split_tags=split_tags, split_plevs=split_plevs)
    tn_split[n + 1, :] .= H_split
    psi = insert_projectors(H_split, psi_split, projs)
    psi_split = split_links(psi; split_tags=split_tags, split_plevs=split_plevs)
  end
  return tn_split, U, Ud
end

struct ContractionAlgorithm{algorithm} end

ContractionAlgorithm(s::AbstractString) = ContractionAlgorithm{Symbol(s)}()

macro ContractionAlgorithm_str(s)
  :(ContractionAlgorithm{$(Expr(:quote, Symbol(s)))})
end

Base.@kwdef struct BoundaryMPS
  top::Vector{MPS}=MPS[]
  bottom::Vector{MPS}=MPS[]
  left::Vector{MPS}=MPS[]
  right::Vector{MPS}=MPS[]
end

struct BoundaryMPSDir{dir} end

BoundaryMPSDir(s::AbstractString) = BoundaryMPSDir{Symbol(s)}()

macro BoundaryMPSDir_str(s)
  :(BoundaryMPSDir{$(Expr(:quote, Symbol(s)))})
end

function contract_approx(tn::Matrix{ITensor}, alg::ContractionAlgorithm"boundary_mps",
                         dirs::BoundaryMPSDir"top_to_bottom"; cutoff, maxdim)
  nrows, ncols = size(tn)
  boundary_mps = Vector{MPS}(undef, nrows - 1)
  x = MPS(tn[1, :])
  boundary_mps[1] = orthogonalize(x, ncols)
  for nrow in 2:(nrows - 1)
    A = MPO(tn[nrow, :])
    x = contract(A, x; cutoff=cutoff, maxdim=maxdim)
    boundary_mps[nrow] = orthogonalize(x, ncols)
  end
  return BoundaryMPS(top=boundary_mps)
end

function contract_approx(tn::Matrix{ITensor}, alg::ContractionAlgorithm"boundary_mps",
                         dirs::BoundaryMPSDir"bottom_to_top"; cutoff, maxdim)
  tn = rot180(tn)
  tn = reverse(tn; dims=2)
  boundary_mps_top = contract_approx(tn, alg, BoundaryMPSDir"top_to_bottom"(); cutoff=cutoff, maxdim=maxdim)
  return BoundaryMPS(bottom=reverse(boundary_mps_top.top))
end

function contract_approx(tn::Matrix{ITensor}, alg::ContractionAlgorithm"boundary_mps",
                         dirs::BoundaryMPSDir"left_to_right"; cutoff, maxdim)
  tn = rotr90(tn)
  boundary_mps = contract_approx(tn, alg, BoundaryMPSDir"top_to_bottom"(); cutoff=cutoff, maxdim=maxdim)
  return BoundaryMPS(left=boundary_mps.top)
end

function contract_approx(tn::Matrix{ITensor}, alg::ContractionAlgorithm"boundary_mps",
                         dirs::BoundaryMPSDir"right_to_left"; cutoff, maxdim)
  tn = rotr90(tn)
  boundary_mps = contract_approx(tn, alg, BoundaryMPSDir"bottom_to_top"(); cutoff=cutoff, maxdim=maxdim)
  return BoundaryMPS(right=boundary_mps.bottom)
end

function contract_approx(tn::Matrix{ITensor}, alg::ContractionAlgorithm"boundary_mps",
                         dirs::BoundaryMPSDir"all"; cutoff, maxdim)
  mps_top = contract_approx(tn, alg, BoundaryMPSDir"top_to_bottom"(); cutoff=cutoff, maxdim=maxdim)
  mps_bottom = contract_approx(tn, alg, BoundaryMPSDir"bottom_to_top"(); cutoff=cutoff, maxdim=maxdim)
  mps_left = contract_approx(tn, alg, BoundaryMPSDir"left_to_right"(); cutoff=cutoff, maxdim=maxdim)
  mps_right = contract_approx(tn, alg, BoundaryMPSDir"right_to_left"(); cutoff=cutoff, maxdim=maxdim)
  return BoundaryMPS(top=mps_top.top, bottom=mps_bottom.bottom, left=mps_left.left, right=mps_right.right)
end

function contract_approx(tn::Matrix{ITensor}, alg::ContractionAlgorithm"boundary_mps";
                         dirs, cutoff, maxdim)
  return contract_approx(tn, alg, BoundaryMPSDir(dirs); cutoff=cutoff, maxdim=maxdim)
end

# Compute the truncation projectors for the network,
# contracting from top to bottom
function contract_approx(tn::Matrix{ITensor}; alg="boundary_mps",
                         dirs="all", cutoff=1e-8, maxdim=maxdim_arg(tn))
  return contract_approx(tn, ContractionAlgorithm(alg); dirs=dirs, cutoff=cutoff, maxdim=maxdim)
end

function insert_projectors(tn::Matrix{ITensor}, boundary_mps::Vector{MPS}; dirs, projector_center)
  return insert_projectors(tn, boundary_mps, BoundaryMPSDir(dirs); projector_center=projector_center)
end

get_itensor(x::MPS, n::Int) = n in 1:length(x) ? x[n] : ITensor()

Base.reverse(x::MPS) = MPS(reverse(x.data))

# From an MPS, create a 1-site projector onto the MPS basis
function projector(x::MPS, projector_center)
  # Gauge the boundary MPS towards the projector_center column
  x = orthogonalize(x, projector_center)

  l = commoninds(get_itensor(x, projector_center - 1), x[projector_center])
  r = commoninds(get_itensor(x, projector_center + 1), x[projector_center])

  uₗ = x[1:(projector_center - 1)]
  uᵣ = reverse(x[(projector_center + 1):end])
  nₗ = length(uₗ)
  nᵣ = length(uᵣ)

  uₗᴴ = dag.(uₗ)
  uᵣᴴ = dag.(uᵣ)

  uₗ′ = reverse(prime.(uₗ))
  uᵣ′ = reverse(prime.(uᵣ))

  if !isempty(uₗ′)
    uₗ′[1] = replaceinds(uₗ′[1], l' => l)
  end
  if !isempty(uᵣ′)
    uᵣ′[1] = replaceinds(uᵣ′[1], r' => r)
  end

  Pₗ = vcat(uₗᴴ, uₗ′)
  Pᵣ = vcat(uᵣᴴ, uᵣ′)
  return Pₗ, Pᵣ
end

default_projector_center(tn::Matrix) = (:, (size(tn, 2) + 1) ÷ 2)

function split_network(tn::Matrix{ITensor}; projector_center=default_projector_center(tn))
  @assert (projector_center[1] == :)
  tn_split = copy(tn)
  nrows, ncols = size(tn)
  @assert length(projector_center) == 2
  @assert projector_center[1] == Colon()
  projector_center_cols = projector_center[2]
  for ncol in 1:ncols
    if ncol ∉ projector_center_cols
      tn_split[:, ncol] = split_links(tn_split[:, ncol])
    end
  end
  return tn_split
end

function insert_projectors(tn::Matrix{ITensor}, boundary_mps::Vector{MPS},
                           dirs::BoundaryMPSDir"top_to_bottom"; projector_center)
  nrows, ncols = size(tn)
  @assert length(boundary_mps) == nrows - 1
  @assert all(x -> length(x) == ncols, boundary_mps)
  @assert length(projector_center) == 2
  @assert (projector_center[1] == :)
  projector_center_cols = projector_center[2]
  projectors = [projector(x, projector_center_cols) for x in boundary_mps]
  projectors_left = first.(projectors)
  projectors_right = last.(projectors)
  tn_split = split_network(tn; projector_center=projector_center)
  return tn_split, projectors_left, projectors_right
end

function insert_projectors(tn, boundary_mps::BoundaryMPS; center, projector_center=nothing)
  top_bottom = true
  if (center[1] == :)
    center = reverse(center)
    tn = rotr90(tn)
    top_bottom = false
  end

  @assert (center[2] == :)
  center_row = center[1]

  if isnothing(projector_center)
    projector_center = default_projector_center(tn)
  end
  @assert (projector_center[1] == :)

  # Approximately contract the tensor network.
  # Outputs a Vector of boundary MPS.
  #boundary_mps_top = contract_approx(tn; alg="boundary_mps", dirs="top_to_bottom", cutoff=cutoff, maxdim=maxdim)
  #boundary_mps_bottom = contract_approx(tn; alg="boundary_mps", dirs="bottom_to_top", cutoff=cutoff, maxdim=maxdim)

  boundary_mps_top = top_bottom ? boundary_mps.top : boundary_mps.left
  boundary_mps_bottom = top_bottom ? boundary_mps.bottom : boundary_mps.right

  # Insert approximate projectors into rows of the network
  boundary_mps_tot = vcat(boundary_mps_top[1:(center_row - 1)], dag.(boundary_mps_bottom[center_row:end]))
  tn_split, projectors_left, projectors_right = insert_projectors(tn, boundary_mps_tot; dirs="top_to_bottom", projector_center=projector_center)

  if !top_bottom
    tn_split = rotl90(tn_split)
  end
  return tn_split, projectors_left, projectors_right
end

# TODO: Delete this. Better to call contract_approx seperately
# and cache the results so they can be reused in different
# calls to `insert_projectors` that have different values of `center`
function insert_projectors(tn; center, projector_center=nothing, maxdim, cutoff)
  # For now, only support contracting from top-to-bottom
  # and bottom-to-top down to a specified row of the network
  top_bottom = true
  if (center[1] == :)
    center = reverse(center)
    tn = rotr90(tn)
    top_bottom = false
  end

  @assert (center[2] == :)
  center_row = center[1]

  if isnothing(projector_center)
    projector_center = (:, (size(tn, 2) + 1) ÷ 2)
  end
  @assert (projector_center[1] == :)

  # Approximately contract the tensor network.
  # Outputs a Vector of boundary MPS.
  boundary_mps_top = contract_approx(tn; alg="boundary_mps", dirs="top_to_bottom", cutoff=cutoff, maxdim=maxdim)
  boundary_mps_bottom = contract_approx(tn; alg="boundary_mps", dirs="bottom_to_top", cutoff=cutoff, maxdim=maxdim)
  # Insert approximate projectors into rows of the network
  boundary_mps = vcat(boundary_mps_top[1:(center_row - 1)], dag.(boundary_mps_bottom[center_row:end]))
  tn_split, projectors_left, projectors_right = insert_projectors(tn, boundary_mps; dirs="top_to_bottom", projector_center=projector_center)

  if !top_bottom
    tn_split = rotl90(tn_split)
  end
  return tn_split, projectors_left, projectors_right
end

function filterneighbors(f, tn, n)
  neighbors_tn = keytype(tn)[]
  tnₙ = tn[n]
  for m in keys(tn)
    if f(n, m) && hascommoninds(tnₙ, tn[m])
      push!(neighbors_tn, m)
    end
  end
  return neighbors_tn
end
neighbors(tn, n) = filterneighbors(≠, tn, n)
inneighbors(tn, n) = filterneighbors(>, tn, n)
outneighbors(tn, n) = filterneighbors(<, tn, n)

function mapinds(f, ::typeof(linkinds), tn)
  tn′ = copy(tn)
  for n in keys(tn)
    for nn in neighbors(tn, n)
      commonindsₙ = commoninds(tn[n], tn[nn])
      tn′[n] = replaceinds(tn′[n], commonindsₙ => f(commonindsₙ))
    end
  end
  return tn′
end

ITensors.prime(::typeof(linkinds), tn, args...) = mapinds(x -> prime(x, args...), linkinds, tn)
ITensors.addtags(::typeof(linkinds), tn, args...) = mapinds(x -> addtags(x, args...), linkinds, tn)

# Return a dictionary from a site to a combiner
function combiners(::typeof(linkinds), tn)
  Cs = Dict(keys(tn) .=> (ITensor[] for _ in keys(tn)))
  for n in keys(tn)
    for nn in inneighbors(tn, n)
      commonindsₙ = commoninds(tn[n], tn[nn])
      C = combiner(commonindsₙ)
      push!(Cs[n], C)
      push!(Cs[nn], dag(C))
    end
  end
  return Cs
end

function contraction_cache_top(tn, boundary_mps::Vector{MPS}, n)
  tn_cache = fill(ITensor(1.0), size(tn))
  for nrow in 1:size(tn, 1)
    if nrow == n
      tn_cache[nrow, :] .= boundary_mps[nrow]
    elseif nrow > n
      tn_cache[nrow, :] = tn[nrow, :]
    end
  end
  return tn_cache
end

function contraction_cache_top(tn, boundary_mps::Vector{MPS})
  tn_cache_top = Vector{typeof(tn)}(undef, length(boundary_mps))
  for n in 1:length(boundary_mps)
    tn_cache_top[n] = contraction_cache_top(tn, boundary_mps, n)
  end
  return tn_cache_top
end

function contraction_cache_bottom(tn, boundary_mps::Vector{MPS})
  tn_top = rot180(tn)
  tn_top = reverse(tn_top; dims=2)
  boundary_mps_top = reverse(boundary_mps)
  tn_cache_top = contraction_cache_top(tn_top, boundary_mps_top)
  tn_cache = reverse.(tn_cache_top; dims=2)
  tn_cache = rot180.(tn_cache)
  return tn_cache
end

function contraction_cache_left(tn, boundary_mps::Vector{MPS})
  tn_top = rotr90(tn)
  tn_cache_top = contraction_cache_top(tn_top, boundary_mps)
  tn_cache = rotl90.(tn_cache_top)
  return tn_cache
end

function contraction_cache_right(tn, boundary_mps::Vector{MPS})
  tn_bottom = rotr90(tn)
  tn_cache_bottom = contraction_cache_bottom(tn_bottom, boundary_mps)
  tn_cache = rotl90.(tn_cache_bottom)
  return tn_cache
end

function contraction_cache(tn, boundary_mps)
  cache_top = contraction_cache_top(tn, boundary_mps.top)
  cache_bottom = contraction_cache_bottom(tn, boundary_mps.bottom)
  cache_left = contraction_cache_left(tn, boundary_mps.left)
  cache_right = contraction_cache_right(tn, boundary_mps.right)
  return (top=cache_top, bottom=cache_bottom, left=cache_left, right=cache_right)
end

function insert_gauge(tn::NamedTuple, gauge)
  return map(x -> insert_gauge.(x, (gauge,)), tn)
end

function insert_gauge(tn, gauge)
  tn′ = copy(tn)
  for n in keys(gauge)
    for g in gauge[n]
      if hascommoninds(tn′[n], g)
        tn′[n] *= g
      end
    end
  end
  return tn′
end

function boundary_mps_top(tn)
  return [MPS(tn[n][n, :]) for n in 1:length(tn)]
end

function boundary_mps_bottom(tn)
  tn_top = rot180.(tn)
  tn_top = reverse.(tn_top; dims=2)
  boundary_mps = boundary_mps_top(tn_top)
  return reverse(boundary_mps)
end

function boundary_mps_left(tn)
  tn_top = rotr90.(tn)
  return boundary_mps_top(tn_top)
end

function boundary_mps_right(tn)
  tn_bottom = rotr90.(tn)
  return boundary_mps_bottom(tn_bottom)
end

function boundary_mps(tn::NamedTuple)
  _boundary_mps_top = boundary_mps_top(tn.top)
  _boundary_mps_bottom = boundary_mps_bottom(tn.bottom)
  _boundary_mps_left = boundary_mps_left(tn.left)
  _boundary_mps_right = boundary_mps_right(tn.right)
  return BoundaryMPS(top=_boundary_mps_top, bottom=_boundary_mps_bottom, left=_boundary_mps_left, right=_boundary_mps_right)
end

# Return a network that when contracted equals the
# squared norm of the input network
function sqnorm(ψ::Matrix{ITensor})
  # Square the tensor network
  ψᴴ = addtags(linkinds, dag.(ψ), "bra")
  ψ′ = addtags(linkinds, ψ, "ket")
  return mapreduce(vec, vcat, (ψᴴ, ψ′))
end

# Return a network that approximates the squared norm of the
# input network
function sqnorm_approx(ψ::Matrix{ITensor}; center, cutoff, maxdim)
  # Square the tensor network
  ψᴴ = addtags(linkinds, dag.(ψ), "bra")
  ψ′ = addtags(linkinds, ψ, "ket")
  # TODO: implement contract(commoninds, ψ′, ψᴴ)
  tn = ψ′ .* ψᴴ

  # Contract in every direction
  combiner_gauge = combiners(linkinds, tn)
  tnᶜ = insert_gauge(tn, combiner_gauge)
  boundary_mpsᶜ = contract_approx(tnᶜ; maxdim=maxdim, cutoff=cutoff)

  tn_cacheᶜ = contraction_cache(tnᶜ, boundary_mpsᶜ)
  tn_cache = insert_gauge(tn_cacheᶜ, combiner_gauge)
  _boundary_mps = boundary_mps(tn_cache)

  #
  # Insert projectors horizontally (to measure e.g. properties
  # in a row of the network)
  #

  tn_projected = insert_projectors(tn, _boundary_mps; center=center)
  tn_split, Pl, Pr = tn_projected

  ψᴴ_split = split_network(ψᴴ)
  ψ′_split = split_network(ψ′)

  Pl_flat = reduce(vcat, Pl)
  Pr_flat = reduce(vcat, Pr)
  return mapreduce(vec, vcat, (ψᴴ_split, ψ′_split, Pl_flat, Pr_flat))
end


text:
# TO RUN:
# include("main.jl")
#
# AUTHOR
# Dr. Umberto Saetti, Assistant Professor, Department of Aerospace Engineering,
# Auburn university
# saetti@auburn.edu
#
# DATE
# 1/3/2022
#
# DESCRIPTION
# Nonlinear simulation model of the rigid-body dynamics of an F-16
# aircraft. The simulation model is based on Refs 1 and 2. The trim,
# linearization routines, and overall architecture of the simulation
# follows the teachings of Dr. Joe Horn at Penn State. When using, please cite
# Ref. 3.
#
# REFERENCES
# 1) Stevens, B. L., Lewis, F. L., and Johnson E. N., “Aircraft Control
#    and Simulation: Dynamics, Controls, Design, and Autonomous Systems,”
#    Wiley, 3rd Edition, 2015.
#    ISBN: 978-1-118-87098-3
# 2) Garza, F. R., and Morelli, E. A., "A Collection of Nonlinear Aircraft
#    Simulations in MATLAB," NASA TM-2003-212145, 2003.
# 3) Saetti, U., and Horn, J. F., "Flight Simulation and Control using the Julia
#    Language", AIAA Scitech Forum, San Diego, CA, Jan 3-7, 2022.
#    DOI: https://arc.aiaa.org/doi/10.2514/6.2022-2354.
#
# STATES             | NINDICES    | DESCRIPTION
# _________________________________________________________________________
# V_T                | 1           | airspeed [ft/s]
# alpha              | 2           | angle of attach [rad]
# beta               | 3           | sideslip [rad]
# phi theta psi      | 4  5  6     | Euler angles [rad]
# p q r              | 7  8  9     | angula rates [rad]
# x                  | 10          | north position [ft]
# y                  | 11          | east position [ft]
# z                  | 12          | altitude [ft]
# P_a                | 13          | power [lb-ft^2/s]
#
# CONTROL INPUTS     | INDICED     | DESCRIPTION     | RANGE
# _________________________________________________________________________
# delta_t            | 1           | throttle [0-1]  |  (0, 1)
# delta_a            | 2           | aileron [deg]   |  (-21.5, 21.5) deg
# delta_e            | 3           | elevator [deg]  |  (-25, 25) deg
# delta_r            | 4           | rudder [deg]    |  (-30, 30) deg
#
# -------------------------------------------------------------------------

# include packages
#using DifferentialEquations
using LinearAlgebra
using SparseArrays
using ControlSystems
using MAT
using Interpolations
using ApproxFun
using Printf
#ENV["MPLBACKEND"]="tkagg"
using Plots
using LaTeXStrings
using JLD
#using PyPlot
#using GR

@printf("\n                     J-F16           \n")
@printf("\n          Author: Dr. Umberto Saetti\n")
@printf("\n-----------------------------------------------\n")

# include functions
include("atan2.jl")
include("linearize.jl")
include("trimmer.jl")
#include("finp_F16_loop.jl")
include("THRUST.jl")
include("DAMP.jl")
include("DNDR.jl")
include("DNDA.jl")
include("DLDR.jl")
include("DLDA.jl")
include("CN.jl")
include("CL.jl")
include("CM.jl")
include("CZ.jl")
include("CY.jl")
include("CX.jl")
include("RTAU.jl")
include("TGEAR.jl")
include("PDOT.jl")
include("ADC.jl")
include("fix.jl")
include("F16.jl")
include("rk4.jl")
include("finp.jl")
include("simulate.jl")
include("simulate_DI.jl")
include("F16_DI.jl")
# load F-16 constants
include("F16_constants.jl")

# ------------------------ TRIM FLIGHT DYNAMICS MODEL --------------------------

# set velocity [ft/s]
speed=600.0
# vertical flight path angle [deg]
gamma=0.0
# horizontal flight path angle [deg]
chi=0.0
# set initial altitude [ft]
alt=2000.0
# set initial North-South and East-West position [ft]
p_north=0.0
p_east=0.0
# turn rate [deg/s]
turnrate=0.0
# pull-up rate [deg/s]
pullup=0.0
# roll rate [deg/s]
rollrate=0.0
# initial guess for state and control vectors
phi_est=turnrate*D2R*speed/32.17
pndot=speed*cos(gamma*D2R)*cos(chi*D2R)
pedot=speed*cos(gamma*D2R)*sin(chi*D2R)
hdot=speed*sin(gamma*D2R)
x0=[speed; 0.0; 0.0; phi_est; 0; chi*D2R; 0.0; 0.0; 0.0; p_north; p_east; alt;
    0.0]
u0=[0.2; 0; 0; 0]
targ_des=[0; 0; 0; rollrate*D2R; pullup*D2R; turnrate*D2R; 0; 0; 0; pndot;
          pedot; hdot; 0; 0]

# trim aircraft
x0, u0, iter = trimmer(F16!,x0,u0,targ_des)
# linearize
A, Bmat, C, D = linearize(F16!,x0,u0)
# partition in longitudinal and lateral flight dynamics
Alon=A[[1,2,5,8],[1,2,5,8]]
Alat=A[[3,4,7,9],[3,4,7,9]]

# ----------------------------- SPECTRAL ANALYSIS ------------------------------

# eigenvalues
eigs=eigvals(A)
eigsLon=eigvals(Alon)
eigsLat=eigvals(Alat)
# store eigenvalues
save("eigs_J-SimpleHel.jld","eigs",eigs)
# plot eigenvalues
gr(size=(800,600))
plt=plot(real(eigs),imag(eigs),seriestype=:scatter,
      label="", legendfontsize=12, legend=:topleft,
      marker = (:circle, :royalblue, 6), markerstrokecolor = :royalblue)
xaxis!("Real",xguidefontsize=14)
yaxis!("Imag",yguidefontsize=14)
display(plt)
savefig(".\\Plots\\eigs_600fts.png")
savefig(".\\Plots\\eigs_600fts.svg")

# ------------------------------ TIME SIMULATION -------------------------------

# time step [s]
dt=0.01
# length of simulation [s]
Tsim=8
# run open-loop simulation
state_OL, time_OL = simulate(F16!,finp,Tsim,dt,x0,u0)
# number of dynamic inverse states
NDISTATES=9
# initial state vector of closed-loop simulation
x0_DI=[x0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; u0[2]; u0[3]; u0[4]]
# run closed-loop simulation
state_CL, time_CL = simulate_DI(F16_DI!,finp,Tsim,dt,x0,u0)

# plot attitude
gr(size=(800,600))
i=4
p1=plot(time_OL,state_OL[i,:]*180/pi,label="Open Loop",legend=:bottomright, legendfontsize=12,line=(3, :royalblue, :solid))
plot!(time_OL,-state_CL[i,:]*180/pi,label="Closed Loop",legend=:bottomright, legendfontsize=12,line=(3, :brown3, :dash))
yaxis!("φ [deg]",yguidefontsize=14)
xlims!((0,8))
i=5
p2=plot(time_OL,state_OL[i,:]*180/pi,label="",line=(3, :royalblue, :solid))
plot!(time_OL,state_CL[i,:]*180/pi,label="",legend=:bottomright, legendfontsize=12,line=(3, :brown3, :dash))
yaxis!("θ [deg]",yguidefontsize=14)
xlims!((0,8))
i=6
p3=plot(time_OL,state_OL[i,:]*180/pi,label="",line=(3, :royalblue, :solid))
plot!(time_OL,-state_CL[i,:]*180/pi,label="",legend=:bottomright, legendfontsize=12,line=(3, :brown3, :dash))
xaxis!("Time [s]",xguidefontsize=14)
yaxis!("ψ [deg]",yguidefontsize=14)
xlims!((0,8))
p4=plot(p1, p2, p3, layout=(3,1), leftmargin=3Plots.mm)
display(p4)
savefig(".\\Plots\\att_F16.svg")

# plot angular rates
gr(size=(800,600))
i=7
p5=plot(time_OL,state_OL[i,:],label="",line=(3, :royalblue, :solid))
plot!(time_OL,-state_CL[i,:],label="",line=(3, :brown3, :dash))
yaxis!("p [rad/s]",yguidefontsize=14)
xlims!((0,8))
i=8
p6=plot(time_OL,state_OL[i,:],label="",line=(3, :royalblue, :solid))
plot!(time_OL,state_CL[i,:],label="",legend=:bottomright, legendfontsize=12,line=(3, :brown3, :dash))
yaxis!("q [rad/s]",yguidefontsize=14)
xlims!((0,8))
i=9
p7=plot(time_OL,state_OL[i,:],label="Open Loop",legend=:bottomright, legendfontsize=12,line=(3, :royalblue, :solid))
plot!(time_OL,-state_CL[i,:],label="Closed Loop",line=(3, :brown3, :dash))
xaxis!("Time [s]",xguidefontsize=14)
yaxis!("r [rad/s]",yguidefontsize=14)
xlims!((0,8))
p8=plot(p5, p6, p7, layout=(3,1), leftmargin=3Plots.mm)
display(p8)
savefig(".\\Plots\\ang_F16.svg")


text:
using ABMExamples
using Test
using Statistics: mean

@testset "ABMExamples.jl" begin
    @testset "SchellingsSegregation.jl" begin
        schelling_data, schelling_filename = run_schelling_model!(20,"schelling")
        @show mean(schelling_data.sum_mood)
        @test mean(schelling_data.sum_mood) == 274.0
    end
    @testset "ForestFire.jl" begin
        _,ff_data = run_forest_fire_model!("forest_new.mp4")
        rounded = round(mean(ff_data.burnt_percentage), digits=2)
        @test rounded == 0.36
    end
    @testset "HKOpinionDynamics.jl" begin
        all_data = hk_model_run_and_plot!()
        # we are testing only ϵ = 0.3, which is the last in the list
        data_epsilon_03 = all_data[end]
        rounded = round(
            mean(data_epsilon_03.new_opinion),
            digits=2
        )
        @test rounded == 0.54
    end
    @testset "Flocking.jl" begin
        run_flocking_example!()
        @test isfile("flocking.mp4") == true
    end
    @testset "PredatorPrey.jl" begin
        _,adf,mdf = run_predatorprey_model!()
        rounded_sheep = round(mean(adf.count_sheep), digits=2)
        rounded_wolves = round(mean(adf.count_wolves), digits=2)
        @test rounded_sheep == 62.99
        @test rounded_wolves == 15.09
    end
end


text:
module AMDGPUnative

using LLVM, LLVM.Interop
using InteractiveUtils
using HSARuntime
using Adapt
using Libdl

const configured = HSARuntime.configured

# where the ROCm-Device-Libs bitcode goes
include(joinpath(@__DIR__, "..", "deps", "deps.jl"))
const device_libs_path = joinpath(@__DIR__, "..", "deps", "usr", "lib")

# needs to be loaded _before_ the compiler infrastructure, because of generated functions
include(joinpath("device", "tools.jl"))
include(joinpath("device", "pointer.jl"))
include(joinpath("device", "array.jl"))
include(joinpath("device", "gcn_intrinsics.jl"))
include(joinpath("device", "runtime_intrinsics.jl"))

include("execution_utils.jl")
include("compiler.jl")
include("execution.jl")

include("reflection.jl")

function __init__()
    check_deps()
end

end # module


text:
export DMRGOpts

"""
Options for the DMRG minimization algorithm.

Fields:
    algorithm::String: DMRG algorithm used for minimizing the state. Can be:
        "DMRG1" for one site DMRG, "DMRG2" for two site DMRG, and "DMRG3S" for
        strictly single site DMRG with subspace expansion.
    nsweeps::Int: number of sweeps.
    maxm::Vector{Int}: maximum bond dimension at every sweep.
    cutoff::Vector{Float64}: cutoff in SVD factorization at every sweep.
    dmrg_tol::Float64: stop the algorithm when the change in the variance of the
        state is less than dmrg_tol.
    show_trace::Int: output information at every step of the minmization:
        - 0: no info given.
        - 1: energy, variance and their variation with respect to their last
            value after every right + left sweep.
        - 2: energy and size of the local Hamiltonian at every step of every
            sweep.
    lanczos_iters::Int: number of maximum iterations in Arpack solver for
        minimizing local Hamiltonians.
    α::Float64: noise term in DMRG3S.
"""
struct DMRGOpts
    algorithm::String
    nsweeps::Int
    maxm::Vector{Int}
    cutoff::Vector{Float64}
    dmrg_tol::Float64
    show_trace::Int
    lanczos_iters::Int
    α::Float64
end

"""
    DMRGOpts(algorithm::String,
             maxm::Vector{Int},
             cutoff::Vector{Float64};
             dmrg_tol::Float64=1e-6,
             show_trace::Int=1,
             lanczos_iters::Int=6,
             α::Float64=1e-6)

Make DMRG options.
"""
function DMRGOpts(algorithm::String,
                  maxm::Vector{Int},
                  cutoff::Vector{Float64};
                  dmrg_tol::Float64=1e-6,
                  show_trace::Int=1,
                  lanczos_iters::Int=6,
                  α::Float64=1e-6)
    length(maxm) != length(cutoff) && throw(
        DimensionMismatch("maxm and cutoff have different dimensions")
    )
    return DMRGOpts(algorithm, length(maxm), maxm, cutoff,
                    dmrg_tol, show_trace, lanczos_iters, α)
end

"""
    DMRGOpts(algorithm::String,
             maxm::Vector{Int},
             cutoff::Float64;
             kwargs...)

Make DMRG options.
"""
function DMRGOpts(algorithm::String,
                  maxm::Vector{Int},
                  cutoff::Float64;
                  kwargs...)
    return DMRGOpts(algorithm, maxm, fill(cutoff, length(maxm)); kwargs...)
end

"""
    DMRGOpts(algorithm::String,
             maxm::Vector{Int};
             kwargs...)

Make DMRG options.
"""
function DMRGOpts(algorithm::String,
                  maxm::Vector{Int};
                  kwargs...)
    cutoff = get(kwargs, :cutoff, 1e-8)
    return DMRGOpts(algorithm, maxm, cutoff; kwargs...)
end

import Base.display

function Base.display(dmrg_opts::DMRGOpts)
    println("DMRG Options")
    println("  Algorithm: $(dmrg_opts.algorithm)")
    println("  Number of sweeps: $(dmrg_opts.nsweeps)")
    println("  DMRG tolerance: $(dmrg_opts.dmrg_tol)")
    println("  Show trace: $(dmrg_opts.show_trace)")
    println("  Lanczos iterations: $(dmrg_opts.lanczos_iters)")
    println("  maxm   cutoff")
    for i in eachindex(dmrg_opts.maxm)
        println("  $(dmrg_opts.maxm[i])     $(dmrg_opts.cutoff[i])")
    end
end


text:
abstract type AbstractPitotTube end

struct Pitot <: AbstractPitotTube
    "Type of Pitot tube (Prandtl, Cole, etc)"
    type::String
    "Tag of the tube - book keeping"
    tag::String
    "Discharge Coefficient"
    Cd::Float64
end
Base.broadcastable(P::Pitot) = Ref(P) # It doesn't broadcast

raw"""
    `Pitot(;Cd=1.0, type="Prandtl", tag="")`

Standard interface for single degree of freedom Pitot tubes.

The velocity of a Pitot tube is given by

``
V = C_d ⋅ \sqrt{2ΔP/ρ}
``
"""
Pitot(;Cd=1.0, type="Prandtl", tag="") = Pitot(type, tag, Cd)

"Return the discharge coefficient of a Pitot tube"
dischcoef(p::Pitot) = p.Cd

"Calculates the velocity measured by a Pitot tube"
velocity(ΔP, ρ=1.09) = sqrt(2ΔP/ρ)
velocity(p::Pitot, ΔP, ρ=1.09) = p.Cd * sqrt(2ΔP / ρ)

(P::Pitot)(ΔP, ρ=1.09) = velocity(P, ΔP, ρ)

"Calculate the dynamic pressure measured by the Pitot tube"
dynamicpress(ΔP) = ΔP
dynamicpress(p::Pitot, ΔP) = p.Cd^2 * ΔP

"Return the pressure differential required by a Pitot tube to achieve a given velocity"
pressdiff(V, ρ=1.09) = ρ*V^2 / 2
pressdiff(p::Pitot, V, ρ=1.09) = ρ*V^2/(2*p.Cd^2)


text:
function ldl_decomp(A)
    # assumption: the matrix A is positive definite (Hermitian)
    # not explicitly checked but assumed: n_row == n_col
    # only the lower left part is used from A

    (n_row, n_col) = size(A)
    L = UnitLowerTriangular( zeros(n_row, n_col) )
    d = zeros(n_col)

    # initial d element
    @inbounds d[1] = A[1,1]
    @inbounds for jj = 1:(n_col-1)

        # determine the next column of L
        @inbounds for ii = (jj+1):n_row
            internal_sum = zero(eltype(A))
            @inbounds @simd for kk = 1:(jj-1)
                @fastmath internal_sum += d[kk] * L[jj,kk] * L[ii,kk]
            end
            @fastmath L[ii, jj] = (A[ii, jj] - internal_sum)/ d[jj]
        end

        # determine the next d element
        diag_sum = zero(eltype(A))
        @inbounds @simd for kk = 1:jj
            @fastmath diag_sum += d[kk] * L[jj+1, kk]^2
        end
        @fastmath d[jj+1] = A[jj+1, jj+1] - diag_sum
    end

    return (L, d)
end

function linsolve_ldl(A, b)

    (L, d) = ldl_decomp(A)
    y = triangular_solve_lower_colwise(L, b)
    @fastmath y ./= d
    x = triangular_solve_upper_colwise_transposed(L, y)

    return x
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
"""
Compute the envelope of a signal along its centroid frequency (see STFT.centroidfrequency)
Parameters:
    -- data: The signal to find the STFT of
    -- windowFunc: An array sampled from the window function to be applied before taking the FFTs. The length of this function defines the window length
    -- windowStep: The step ( in samples) to move the window between FFTs
    -- low: The lowest harmonic to include in the output
    -- high: The highest harmonic to include in the output
    -- transformFunction: A function to apply to the windowed, FFT'd pieces of _data_
"""
function envelope_along_centroid{dataType,windowType}(data::AbstractArray{dataType, 1}, windowFunc::AbstractArray{windowType, 1}, windowStep::Integer, low::Integer, high::Integer, transformFunction::Function, bandWidth::Integer)
    halfBandWidth = Int(round(bandWidth/2))
    env_along_path = stft_processed(data, windowFunc, windowStep, low, high, transformFunction, (dat, i)->(cent=Int(round(centroid(dat)));mean(view(dat, max(cent-halfBandWidth, 1):min(cent+halfBandWidth, length(dat))))), Float64)
    return env_along_path
end

function envelope_along_centroid{dataType}(data::AbstractArray{dataType, 2}, bandWidth::Integer)
    envLength = size(data, 2)
    envelope = Array{dataType, 1}(envLength)

    for timeSlice in 1:envLength
	centroidIndex = Int(round(centroid(data[:, timeSlice])))
	envelope[timeSlice] = mean(data[(centroidIndex-bandWidth):(centroidIndex+bandWidth), i])
    end

    return envelope
end

function envelope{T<:Any, idxType<:Integer}(data::AbstractArray{T,2}, path::AbstractArray{idxType, 1}, bandWidth::Integer)
  nRows, nCols = size(data)
  env = Array{T, 1}(nCols)

  for colNum in 1:nCols
    env[colNum] = mean(view(data, max(path[colNum]-Int(floor(bandWidth/2)), 1):min(path[colNum]+Int(ceil(bandWidth/2)), nRows), colNum))
  end

  return env
end

export envelope
export envelope_along_centroid


text:
# load libraries and code
using PyPlot
include("schelling_model.jl")

# initialize world
world = initialize_world(100, 0.7)

# show initial state
PyPlot.clf()
PyPlot.subplot(1, 2, 1)
PyPlot.imshow(world)

# execute simulation
new_world = simulate(world, 1000, 0.9)

# show state after simulation
PyPlot.subplot(1, 2, 2)
PyPlot.imshow(new_world)
display(gcf())

max_stars_repo_name: JohannesNakayama/schelling-model

max_issues_repo_name: JohannesNakayama/schelling-model

max_forks_repo_name: JohannesNakayama/schelling-model


text:
"""
Various utility functions for using across models, controllers and views.
"""
module Helpers

using Genie, Sessions, Router, URIParser, Logger, HttpServer, WebSockets

export session, request, response, flash, number_of_pages, paginated_uri, var_dump, wsclient, flash_has_message

"""
    session() :: Sessions.Session
    session(params::Dict{Symbol,Any}) :: Sessions.Session

Returns the `Session` object associated with the current HTTP request.
"""
function session() :: Sessions.Session
  session(Router._params_())
end
function session(params::Dict{Symbol,Any}) :: Sessions.Session
  if haskey(params, Genie.PARAMS_SESSION_KEY)
    return params[Genie.PARAMS_SESSION_KEY]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_SESSION_KEY) key"
    Logger.log(msg, :err)
    error(msg)
  end
end

"""
    request() :: HttpServer.Request
    request(params::Dict{Symbol,Any}) :: HttpServer.Request

Returns the `Request` object associated with the current HTTP request.
"""
function request() :: HttpServer.Request
  request(Router._params_())
end
function request(params::Dict{Symbol,Any}) :: HttpServer.Request
  if haskey(params, Genie.PARAMS_REQUEST_KEY)
    return params[Genie.PARAMS_REQUEST_KEY]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_REQUEST_KEY) key"
    Logger.log(msg, :err)
    error(msg)
  end
end

"""
    response() :: HttpServer.Response
    response(params::Dict{Symbol,Any}) :: HttpServer.Response

Returns the `Response` object associated with the current HTTP request.
"""
function response() :: HttpServer.Response
  response(Router._params_())
end
function response(params::Dict{Symbol,Any}) :: HttpServer.Response
  if haskey(params, Genie.PARAMS_RESPONSE_KEY)
    return params[Genie.PARAMS_RESPONSE_KEY]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_RESPONSE_KEY) key"
    Logger.log(msg, :err)
    error(msg)
  end
end

"""
    flash(params::Dict{Symbol,Any})

Returns the `flash` dict object associated with the current HTTP request.
"""
function flash()
  flash(Router._params_())
end
function flash(params::Dict{Symbol,Any})
  if haskey(params, Genie.PARAMS_FLASH_KEY)
    return params[Genie.PARAMS_FLASH_KEY]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_FLASH_KEY) key"
    Logger.log(msg, :err)
    error(msg)
  end
end

"""
    flash(value::Any) :: Void
    flash(value::Any, params::Dict{Symbol,Any}) :: Void

Stores `value` on the `flash`.
"""
function flash(value::Any) :: Void
  flash(value, Router._params_())
end
function flash(value::Any, params::Dict{Symbol,Any}) :: Void
  Sessions.set!(session(params), Genie.PARAMS_FLASH_KEY, value)
  params[Genie.PARAMS_FLASH_KEY] = value

  nothing
end

"""

"""
function flash_has_message() :: Bool
  ! isempty(flash())
end

"""
    wsclient(params::Dict{Symbol,Any}) :: WebSockets.WebSocket

Returns the `WebSocket` object associated with the current WS request.
"""
function wsclient(params::Dict{Symbol,Any}) :: WebSockets.WebSocket
  if haskey(params, Genie.PARAMS_WS_CLIENT)
    return params[Genie.PARAMS_WS_CLIENT]
  else
    msg = "Invalid params Dict -- must have $(Genie.PARAMS_WS_CLIENT) key"
    Logger.log(msg, :err)
    error(msg)
  end
end

end


text:
"""
`sgf(yr,r,g)`: State Government Finance data.
```math
\\alpha_{yr,r,g}^{sgf} = \\dfrac{\\bar{sgf}_{yr,r,g}}{\\sum_{r'} \\bar{sgf}_{yr,r',g}}
```
Note: D.C. is not included in the original data set.
Assume its SGF data is equal to Maryland's.
"""
function share_sgf!(d::Dict)
    println("  Calculating sgf(yr,r,g), regional shares of State Government Finance data")
    d[:sgf] = d[:sgf] / transform_over(d[:sgf], :r)
    verify_over(d[:sgf], :r) !== true && @error("SGF shares don't sum to 1.")
    return dropnan!(d[:sgf])
end

max_stars_repo_name: NREL/SLiDE

max_issues_repo_name: NREL/SLiDE

max_forks_repo_name: NREL/SLiDE


text:
function total_fuel_cost(experiment_time::String, trial::Int64, day_range, sys::PowerSystems.System)
    result = collect_ed_results(experiment_time, trial, day_range, 24)

    variable_cost = 0.0
    fixed_cost = 0.0
    for n in names(result["cg"])
        variable_cost += sum(result["cg"][!,n])
        g = PowerSystems.get_component(PowerSystems.ThermalStandard, sys, String(n))
        g_fixed_cost = PowerSystems.get_op_cost(g) |> PowerSystems.get_fixed
        fixed_cost += sum(result["ug"][!,n]*g_fixed_cost)
    end

    total_cost = variable_cost + fixed_cost

    return total_cost
end

function total_load_not_supplied(experiment_time::String, trial::Int64, day_range, sys::PowerSystems.System)
    res = collect_ed_results(experiment_time, trial, day_range, 24)
    return aggregate(res["il_ub"] .- res["pl"], sum)[!,1][1]
end

function total_curtailment(experiment_time::String, trial::Int64, day_range, sys::PowerSystems.System)
    res = collect_ed_results(experiment_time, trial,  day_range, 24)
    return aggregate(res["pw_ub"] .- res["pw"], sum)[!,1][1]
end

function total_computing_time(experiment_time::String, trial::Int64, day_range, model::String)
    res = collect_uc_results(experiment_time, model, trial::Int64, day_range)
    return res["optimizer_log"][!,:solve_time]
end