text:
import LightGraphs
import Polyhedra
import GLPK
import Plots
import Random
import Statistics
import Triangulate

"""
    gen_obstacle(max_side_len)

Creates an obstacle from a given set of points in the unit square.
"""
function gen_obstacle(max_side_len::Real)
    # TODO: Make this function accept a list of points to take the convex hull
    # Faces must be passed in counter-clockwise order
    # Use rational points (convert if necessary). Rid redundant points and compute HRep.
    # Assure obstacles have at least 4 extreme points for pairwise IB to work

    # Julia Notes
    # # "varargs" or slurping up arguments
    # function myAdd(x, y, z...) # the ... here means collect additional arguments into z
    #     if isempty(z)
    #         return x + y
    #     end
    #     return x + myAdd(y, z...)
    #     # note the function calls it self (recursion)
    #     # z is a container; the ... is necessary so that its contents can be passed as arguments to myAdd
    # end

    # This might also be helpful
    # v = Polyhedra.convexhull(map(i -> collect(points[i]), face)...)

    # Remove redundant points with

    point = rand(2)

    v = Polyhedra.convexhull(point, map(x -> x > 1 ? 1 : x, point + max_side_len * rand(2)), map(x -> x > 1 ? 1 : x, point + max_side_len * rand(2)))

    if Polyhedra.npoints(poly) < 4
        # something
    end

    poly = Polyhedra.polyhedron(v, Polyhedra.DefaultLibrary{Rational{Int64}}(GLPK.Optimizer))
    Polyhedra.removevredundancy!(poly)
    Polyhedra.hrep(poly)
    return poly

end

"""
    gen_obstacle_four(point_per_dim)

Creates an obstacle with four extreme points from a grid of points in the unit square.
"""
function gen_obstacle_four(points_per_dim::Int)
    points_per_dim = max(points_per_dim, 6)  # Need at least 6 interior points
    nums = LinRange(0//1, 1//1, points_per_dim+2)[2:end-1]  # Do not want obstacles touching the unit square boundaries

    # Select the bottom left point at random and select the other three points with respect to it
    x_index = rand(2:(length(nums)-4))
    y_index = rand(2:(length(nums)-4))
    bottom_left = [ nums[x_index] ; nums[y_index] ]
    bottom_right = [ nums[x_index + 3 + rand(-1:1)] ; nums[y_index + rand(-1:1)] ]
    top_left = [ nums[x_index + rand(-1:1)] ; nums[y_index + 3 + rand(-1:1)] ]
    top_right = [ nums[x_index + 3 + rand(-1:1)] ; nums[y_index + 3 + rand(-1:1)] ]

    v = Polyhedra.convexhull(top_left, top_right, bottom_left, bottom_right)
    poly = Polyhedra.polyhedron(v, Polyhedra.DefaultLibrary{Rational{Int64}}(GLPK.Optimizer))
    Polyhedra.removevredundancy!(poly)
    # Loop below assures the polyhedron has 4 extreme points (3 can occur under certain combinations) for pairwise IB purposes
    while Polyhedra.npoints(poly) < 4
        top_right = [ nums[x_index + 3 + rand(-1:1)] ; nums[y_index + 3 + rand(-1:1)] ]
        v = Polyhedra.convexhull(top_left, top_right, bottom_left, bottom_right)
        poly = Polyhedra.polyhedron(v, Polyhedra.DefaultLibrary{Rational{Int64}}(GLPK.Optimizer))
        Polyhedra.removevredundancy!(poly)
    end

    Polyhedra.hrep(poly)
    return poly
end

"""
    gen_field(num_obstacles; custom, seed)

Creates an obstacle course, an array of obstacles in the unit square. If
obstacles overlap, the convex hull of their combined extreme points will be taken
instead.
"""
function gen_field(num_obstacles::Int; custom::Bool=false, points_it::Set{Vector{T}}=Set{Vector{Rational}}(), seed::Int=1) where {T}
    if custom
        obstacles = []
        for points in points_it
            push!(obstacles, gen_obstacle(points))
        end
    else
        Random.seed!(seed)
        points_per_dim = 20
        obstacles = map(_n -> gen_obstacle_four(points_per_dim), 1:num_obstacles)
    end

    return remove_overlaps(obstacles)
end

"""
    remove_overlaps(obstacles)

Given an obstacle course, check for overlapping obstacles and replace offending
obstacles with the convex hull of their combined extreme points.
"""
function remove_overlaps(obstacles)
    if !has_overlaps(obstacles)
        return obstacles
    end

    res = []
    fused = []

    for i = 1:length(obstacles)
        if i in fused
            continue
        end

        accumulator = obstacles[i]
        for j = 1:length(obstacles)
            if i != j && !(j in fused)
                other = obstacles[j]

                overlap = Polyhedra.intersect(accumulator, other)
                if !isempty(overlap)
                    accumulator = Polyhedra.convexhull(accumulator, other)
                    Polyhedra.removevredundancy!(accumulator)
                    Polyhedra.removehredundancy!(accumulator)
                    push!(fused, j)
                end
            end
        end

        push!(res, accumulator)
    end

    return remove_overlaps(res)
end

"""
    has_overlaps(obstacles)

Returns true if any obstacles in the set are overlapping.
"""
function has_overlaps(obstacles)::Bool
    for i = 1:length(obstacles)
        for j = 1:length(obstacles)
            if i != j
                if !isempty(Polyhedra.intersect(obstacles[i], obstacles[j]))
                    return true
                end
            end
        end
    end

    return false
end

"""
    find_intersections(obstacles)

Analyzes the halfspaces that make up each obstacle and unit square border, and finds
where they intersect. Returns an array of all intersections, as well as an array
mapping each halfspace to the points of intersection that lie on it, sorted. Points
in the interior of obstacles are stored at the end of the array.
"""
function find_intersections(obstacles)
    res = []

    halfspaces = @pipe map(obstacle -> Polyhedra.hrep(obstacle).halfspaces, obstacles) |> Iterators.flatten(_) |> collect(_)
    unique!(halfspaces)

    boundaries = [
        Polyhedra.HalfSpace([-1//1, 0//1], 0//1),
        Polyhedra.HalfSpace([1//1, 0//1], 1//1),
        Polyhedra.HalfSpace([0//1, -1//1], 0//1),
        Polyhedra.HalfSpace([0//1, 1//1], 1//1)
    ]
    filter!(b -> !(b in halfspaces), boundaries)

    lines = [halfspaces; boundaries]

    # Save indices of horizontal and vertical halfspaces to keep track of which points need
    # the same y- or x- coordinates (necessary to sort correctly)
    vert_ind = []
    hor_ind = []
    for i = 1:length(lines)
        if abs(lines[i].a'[1]) == 1//1 && abs(lines[i].a'[2]) == 0//1 # set tolerance?
            push!(vert_ind, i)
        end
        if abs(lines[i].a'[1]) == 0//1 && abs(lines[i].a'[2]) == 1//1
            push!(hor_ind, i)
        end
    end

    points = fill([], length(lines))    # For storing the points on each halfspace
    for i = 1:length(lines)
        points[i] = []
    end

    # Function for determining if a point should be considered a duplicate of another. If so, the index
    # of the existing one is returned
    function point_duplicates(point, res)
        tol = Rational(1.e-13)
        for index in 1:length(res)
            if abs(point.first - res[index].first) < tol && abs(point.second - res[index].second) < tol
                return true, index
            end
        end
        return false, -1
    end

    tol = Rational(1.e-13)
    x_vert = Dict([]) # keys are the indices in "lines"
    y_hor = Dict([])
    for i = 1:length(lines)
        for j = (i + 1):length(lines)
            A = vcat(lines[i].a', lines[j].a')
            b = [lines[i].β; lines[j].β]
            try
                x = A \ b

                # Assures boundary points do not have floating point errors
                if abs(x[1] - 0//1) < tol
                    x[1] = 0//1
                elseif abs(x[1] - 1//1) < tol
                    x[1] = 1//1
                end

                if abs(x[2] - 0//1) < tol
                    x[2] = 0//1
                elseif abs(x[2] - 1//1) < tol
                    x[2] = 1//1
                end

                if (x[1] >= 0//1 && x[1] <= 1//1 && x[2] >= 0//1 && x[2] <= 1//1)
                    # Assures points on vertical or horizontal halfspaces have the same x- or y- coordinates, resp.
                    if i in vert_ind
                        if !(i in keys(x_vert))
                            x_vert[i] = x[1]
                        else
                            x[1] = x_vert[i]
                        end
                    elseif i in hor_ind
                        if !(i in keys(y_hor))
                            y_hor[i] = x[2]
                        else
                            x[2] = y_hor[i]
                        end
                    end
                    if j in vert_ind
                        if !(j in keys(x_vert))
                            x_vert[j] = x[1]
                        else
                            x[1] = x_vert[j]
                        end
                    elseif j in hor_ind
                        if !(j in keys(y_hor))
                            y_hor[j] = x[2]
                        else
                            x[2] = y_hor[j]
                        end
                    end

                    point = x[1] => x[2]

                    duplicate, index = point_duplicates(point, res)
                    if !(duplicate)
                        # If point is not a duplicate, push it into res and the halfspaces it is in
                        push!(res, point)
                        push!(points[i], point)
                        push!(points[j], point)
                    else
                        # If point is a duplicate, push it into the halfspaces it's in if not added yet
                        if !(res[index] in points[i])
                            push!(points[i], res[index])
                        end
                        if !(res[index] in points[j])
                            push!(points[j], res[index])
                        end
                    end
                end
            catch e
                # No intersections found
            end
        end
    end

    # Find vertices in the interior of an obstacle
    points_vec = res
    inside_indices = filter(i -> begin
        inside = false
        for ob in obstacles
            if Polyhedra.ininterior([points_vec[i].first ; points_vec[i].second], ob)
                inside = true
                # println("Inside Point: $(points_vec[i].first), $(points_vec[i].second)")
                break
            end
        end
        return inside
    end, 1:length(points_vec))

    # Move inside vertices to the end of the points vector, so that they may be removed from the points
    # vector and the lightgraph graph without altering the vertex labels and edges remaining.
    # Removal does not occur here.
    for i = 1:(length(inside_indices))
        temp = points_vec[inside_indices[i]]
        points_vec[inside_indices[i]] = points_vec[end - (i-1)]
        points_vec[end - (i-1)] = temp
    end

    sorted_points = map(list -> sort(list), points)

    # Third return value is for knowing how many inside points there are
    return (points_vec, sorted_points, length(inside_indices))
end

"""
    find_points(obs)
Compute the points corresponding to obstacle extreme points and the four corners
of the unit square.
"""
function find_points(obs)
    obstacle_faces = []   # Store obstacle face vertices (should be ordered)
    points = []   # Store points as Vector of Pairs
    vertex_num = 0
    for ob in obs
        vertices_in_ob = []
        iter = Polyhedra.points(ob)
        # Store x and y coordinates in points
        for point in iter
            push!(points, point[1] => point[2])
            vertex_num += 1
            push!(vertices_in_ob, vertex_num)
        end

        push!(obstacle_faces, vertices_in_ob)
    end

    points = [points; 0//1 => 0//1 ; 0//1 => 1//1; 1//1 => 1//1; 1//1 => 0//1]
    # unique!(points)   # Code above relies on the points being unique

    return points
end

"""
    construct_graph(obs)

Given a list of obstacles with no overlaps, finds the intersections of the lines
that make up each obstacle's halfspaces and constructs a graph with the
intersections as nodes and the halfspaces as edges. Returns a tuple containing
the obstacles, a list of the location of the nodes' coordinates, the graph, a set of
every face of obstacle space, and a set of every face of free space. (Each face is
a clockwise-ordered vector of vertices).
"""
function construct_graph(obs)
    points, mapped, inside_quant = ClutteredEnvPathOpt.find_intersections(obs)

    # Create map from point to its neighbors
    neighbors = ClutteredEnvPathOpt._find_neighbors(points, mapped)
    # Create lightgraph
    graph = ClutteredEnvPathOpt._gen_graph(neighbors)

    # TODO: Tolerances issues due to approximating an irrational with a rational
    # angles[i,j] is the angle point j makes with point i
    angles = fill(Rational(Inf), length(points), length(points))
    for i in 1:length(points)
        for j in 1:length(points)
            angles[i, j] = Rational(atan(points[j].second - points[i].second, points[j].first - points[i].first))
            if angles[i, j] < 0//1
                angles[i, j] += 2 * pi # Rational(2*3.141592653589793)
            end
        end
    end

    function greatest_angle_neighbor(source, last_angle, visited)
        #tol = Rational(1.e-15)
        unvisited = filter(vertex -> !(vertex in visited), neighbors[source])
        neighbor_angles = map(neighbor -> angles[source, neighbor], unvisited)
        # For sorting angles
        for i in 1:length(neighbor_angles)
            if neighbor_angles[i] >= Rational(3.141592653589793)  # Use pi?
                neighbor_angles[i] -= Rational(2*3.141592653589793)
            end
        end
        zipped = @pipe zip(unvisited, neighbor_angles) |> collect(_)
        # zipped is an array of tuples of the form: (unvisited node, angle it makes with source)

        # Two cases to assure we are constructing the most compact face (ie, no face with subfaces)
        if last_angle >= 0//1 && last_angle < Rational(3.141592653589793)
            greater = filter(tup -> (tup[2] > last_angle) && (tup[2] < last_angle + Rational(3.141592653589793)), zipped)
        else
            greater = filter(tup -> (tup[2] > last_angle - Rational(2*3.141592653589793)) && (tup[2] < last_angle - Rational(3.141592653589793)), zipped)
        end

        if isempty(greater)
            return (-1, -1) # bad state
        else
            sort!(greater, by=(tup -> tup[2]))
            destination = greater[end][1]
        end

        return (destination, angles[source, destination])
    end

    faces = []

    for i in 1:length(points)
        start = i

        for j in 1:length(neighbors[start])
            current = neighbors[start][j]
            last_angle = angles[start, current]

            face = [start, current]

            while start != current

                if (length(face) > 2 && (start in neighbors[current]))
                    push!(faces, copy(face))
                    break
                end

                current, last_angle = greatest_angle_neighbor(current, last_angle, face)

                if current == -1
                    break
                end

                push!(face, current)
            end
        end
    end

    # Find irredundant faces
    face_sets = []
    unique_faces = filter(face -> begin
        face_set = Set(face)
        included = false

        if !(face_set in face_sets)
            push!(face_sets, face_set)
            included = true
        end

        return included
    end, faces)

    # Find free faces (by determining which do not overlap obstacles)
    free_faces = filter(face -> begin
        face_set = Set(face)   # face is already a set, luckily taking a Set of a Set doesn't do anything
        # TODO: It's not even used! Comment it and test, then remove

        face_set_v = Polyhedra.convexhull(map(i -> collect(points[i]), face)...)
        face_set_poly = Polyhedra.polyhedron(face_set_v, Polyhedra.DefaultLibrary{Rational{Int64}}(GLPK.Optimizer))
        face_set_overlaps_obs_faces = false
        included = false

        # If the intersection of face and an obstacle is the face itself, than that face is a subface of an obstacle
        tol =  Rational(1.e-10)
        for n in 1:length(obs)  # need to loop in this manner for 'continue' to work properly
            indicator = 0
            intersect_poly = Polyhedra.intersect(face_set_poly, obs[n])
            Polyhedra.vrep(intersect_poly)  # compute V-representation

            if isempty(Polyhedra.points(intersect_poly)) == true
                continue  # intersection is empty
            end

            ext_p_intersect_poly = sort(collect(intersect_poly.vrep.points.points))
            ext_p_face_set_poly = sort(collect(face_set_poly.vrep.points.points))
            if length(ext_p_intersect_poly) != length(ext_p_face_set_poly)
                continue  # intersection is not the face
            end

            for i in 1:length(ext_p_face_set_poly)
                for j in 1:length(ext_p_face_set_poly[i])  # should be 1:2 (x and y coordinates)
                    if abs(ext_p_intersect_poly[i][j] - ext_p_face_set_poly[i][j]) > tol
                        indicator += 1  # the extreme points do not coincide at the x or y coordinate
                        break
                    end
                end
                if indicator != 0
                    break
                end
            end
            if indicator == 0
                face_set_overlaps_obs_faces = true
                break
            end
        end

        if !(face_set_overlaps_obs_faces)
            included = true
        end

        return included
    end, unique_faces)

    # Swap obstacle (sub)faces with whole faces, then remove inside nodes and inside edges
    # CHAD SUBSET APPROACH

    # Remove inside vertices
    if inside_quant > 0
        length_p = length(points)
        for i = 1:inside_quant
            LightGraphs.rem_vertex!(graph, length_p - (i-1))
            pop!(points)
        end
    end

    # Determine the vertices that make up each obstacle's extreme points
    obs_euclid_points = map(ob -> Polyhedra.points(ob), obs)
    obstacle_faces = map(ob_points -> begin
        tol = Rational(1e-13)
        face_points = []
        for ob_point in ob_points
            vertex_index = findfirst(p -> abs(ob_point[1] - p.first) < tol && abs(ob_point[2] - p.second) < tol, points)
            if vertex_index !== nothing
                push!(face_points, vertex_index)
            end
        end
        return face_points
    end, obs_euclid_points)

    # Find all vertices associated with an obstacle face and remove inside edges

    boundary_edges = fill([], length(obs))  # store boundary edges of each obstacle
    for i = 1:length(obs)
        boundary_edges[i] = []
    end

    edges_to_remove = []  # edges crossing through an obstacle are to be removed
    for edge in LightGraphs.edges(graph)
        on_boundary = true
        for i = 1:length(obs)
            if Polyhedra.in([points[edge.src].first; points[edge.src].second], obs[i]) && Polyhedra.in([points[edge.dst].first; points[edge.dst].second], obs[i])
                # Edge is contained in obstacle i
                # println("Found Edge: $(edge.src) => $(edge.dst) with i = $i")
                push!(obstacle_faces[i], edge.src)
                push!(obstacle_faces[i], edge.dst)
                unique!(obstacle_faces[i])
                if Polyhedra.ininterior([(points[edge.src].first + points[edge.dst].first)//2 ; (points[edge.src].second + points[edge.dst].second)//2], obs[i])
                    # Edge crosses through obstacle i
                    #println("Removing edge $(edge.src) => $(edge.dst)")
                    push!(edges_to_remove, (edge.src, edge.dst))
                    on_boundary = false
                end
                if on_boundary   # this can be change to an if-else with the above
                    # println("Adding edge $(edge.src) => $(edge.dst) with i = $i")
                    push!(boundary_edges[i], edge)
                end
                break
            end
        end
    end

    for (a,b) in edges_to_remove
        # println("Removing edge $a => $b")
        LightGraphs.rem_edge!(graph, a, b)
    end

    # for i in 1:length(boundary_edges)
    #     println("Obstacle $i")
    #     for edge in boundary_edges[i]
    #         println("$edge")
    #     end
    # end

    # Now need to order points on obstacle faces counter-clockwise (they're then returned clockwise)

    # Start with left-most point (if multiple, use bottom-most)
    starting_points = []
    for ob_face in obstacle_faces
        x_coor = map(vertex -> points[vertex].first, ob_face)
        min_x = min(x_coor...)
        candidates = filter(vertex -> points[vertex].first == min_x, ob_face)
        if length(candidates) > 1
            y_coor = map(vertex -> points[vertex].second, candidates)
            min_y = min(y_coor...)
            candidates = filter(vertex -> points[vertex].second == min_y, candidates)
        end
        push!(starting_points, candidates[1])
    end

    # println("Starting Points:")
    # for sp in starting_points
    #     println("$sp")
    # end

    obstacle_faces_ordered = []
    for i = 1:length(starting_points)
        sp = starting_points[i]
        sp_adj = []  #  adjacent nodes (should be exactly 2)
        for edge in boundary_edges[i]
            if sp == edge.src
                push!(sp_adj, edge.dst)
            elseif sp == edge.dst
                push!(sp_adj, edge.src)
            end
        end
        sp_n1 = sp_adj[1]
        sp_n2 = sp_adj[2]

        angle_n1 = angles[sp, sp_n1]
        if angle_n1 >= Rational(3.141592653589793)  # for ordering
            angle_n1 -= Rational(2*3.141592653589793)
        end
        angle_n2 = angles[sp, sp_n2]
        if angle_n2 >= Rational(3.141592653589793)
            angle_n2 -= Rational(2*3.141592653589793)
        end

        # Determine counter-clockwise order of the three points
        n1 = sp
        if angle_n1 > angle_n2
            face = [sp_n1; sp; sp_n2]
            n0 = sp_n1
            n2 = sp_n2
        else
            face = [sp_n2; sp; sp_n1]
            n0 = sp_n2
            n2 = sp_n1
        end

        # From current node n2, determine the next node in the face by checking which other node n2
        # is adjacent to (other than n1). Repeat until the next node is n0.
        while n2 != n0
            next_found = false
            for edge in boundary_edges[i]
                if edge.src == n2 && edge.dst != n1
                    n1 = copy(n2)
                    n2 = edge.dst
                    next_found = true
                    break
                elseif edge.dst == n2 && edge.src != n1
                    n1 = copy(n2)
                    n2 = edge.src
                    next_found = true
                    break
                end
            end
            if !next_found
                println("Infinite while loop. Face may not be complete.")
                break
            end
            if n2 != n0
                push!(face, n2)
            end
        end

        push!(obstacle_faces_ordered, face)
    end

    return (obs, points, graph, Set(map(face -> reverse(face), obstacle_faces_ordered)), Set(map(f -> reverse(f), free_faces)))
end

"""
    construct_graph_delaunay(obs)

Generate the constrained Delaunay triangulation partition of the free space.
Returns a tuple containing the obstacles, a list of the location of the nodes' coordinates,
the graph, a set of every face of obstacle space, and a set of every face of free space.
(Each face is a clockwise-ordered vector of vertices).
"""
function construct_graph_delaunay(obs)

    S = []   # Store segments
    obstacle_faces = []   # Store obstacle face vertices (should be ordered)
    points = []   # Store points as Vector of Pairs
    vertex_num = 0
    for ob in obs
        vertices_in_ob = []
        iter = Polyhedra.points(ob)
        # Store x and y coordinates in points
        for point in iter
            push!(points, point[1] => point[2])
            vertex_num += 1
            push!(vertices_in_ob, vertex_num)
        end
        # Push edges of obstacle boundary (assumes they're in order)
        for i = 1:(length(vertices_in_ob)-1)
            push!(S, [vertices_in_ob[i] ; vertices_in_ob[i+1]])
        end
        push!(S, [vertices_in_ob[end] ; vertices_in_ob[1]])

        push!(obstacle_faces, vertices_in_ob)
    end

    points = [points; 0//1 => 0//1 ; 0//1 => 1//1; 1//1 => 1//1; 1//1 => 0//1]
    # unique!(points)   # Code above relies on the points being unique

    # Add boundary segments
    num_points = length(points)
    for i = (num_points-3):(num_points-1)
        push!(S, [i ; i+1])
    end
    push!(S, [num_points ; num_points-3])

    # Initialize triangulate structure
    triin = Triangulate.TriangulateIO()

    # Store points in matrix
    C = zeros(2,length(points))
    for j = 1:length(points)
        C[:,j] = [points[j].first ; points[j].second]
    end
    triin.pointlist = Matrix{Cdouble}(C)

    # Store segments in matrix
    segments = zeros(Int,2,length(S))
    for j=1:length(S)
        segments[:,j] = [S[j][1]; S[j][2]]
    end
    triin.segmentlist = Matrix{Cint}(segments)

    # Add hole list
    hole_center = zeros(2,length(obstacle_faces))
    for (j,face) in enumerate(obstacle_faces)
        x_locations = map(i -> points[i].first, face)
        y_locations = map(i -> points[i].second, face)
        avg_x = Statistics.mean(x_locations)
        avg_y = Statistics.mean(y_locations)
        hole_center[:,j] = [avg_x; avg_y]
    end
    triin.holelist = Matrix{Cdouble}(hole_center)

    # Computed CDT
    (triout, _) = Triangulate.triangulate("pQ", triin)

    triangles = Matrix{Int64}(triout.trianglelist)   # matrix with 3 rows
    # edges = Matrix{Int64}(triout.edgelist)   # matrix with 2 rows

    # Create graph
    graph = LightGraphs.SimpleGraph(length(points))
    # Add edges
    for j = 1:size(triangles,2)
        LightGraphs.add_edge!(graph, triangles[1,j] => triangles[2,j])
        LightGraphs.add_edge!(graph, triangles[1,j] => triangles[3,j])
        LightGraphs.add_edge!(graph, triangles[2,j] => triangles[3,j])
    end
    # Store free faces as vectors (utilizing holes, all triangles are free faces)
    free_faces = []
    for j = 1:size(triangles,2)
        push!(free_faces, triangles[:,j])
    end

    # TODO: Face merging

    return (obs, points, graph, Set(map(face -> face, obstacle_faces)), Set(map(f -> reverse(f), free_faces)))
end

"""
    _find_neighbors(points, mapped)

Given a list of points and an association of halfspace lines to points that lie
on the line, return an array where res[i] <- the neighbors of node i.
"""
function _find_neighbors(points::Vector{Any}, mapped::Vector{Vector{Any}})::Vector{Vector{Any}}
    # This function should only be used in construct_graph(), before "inside" points are removed
    neighbors = Dict()

    for point in points
        neighbors[point] = []

        for colinear in mapped
            if in(point, colinear)
                index = findfirst(p -> p == point, colinear)
                # Do not push the same point more than once
                if index - 1 > 0 && !(colinear[index - 1] in neighbors[point])
                    push!(neighbors[point], colinear[index - 1])
                end

                if index + 1 <= length(colinear) && !(colinear[index + 1] in neighbors[point])
                    push!(neighbors[point], colinear[index + 1])
                end
            end
        end
    end

    # Res will contain the neighbors by their node label
    res = fill([], length(points))
    for i in 1:length(res)
        res[i] = map(point -> findfirst(p -> p == point, points), neighbors[points[i]])
    end

    return res
end

"""
    _gen_graph(neighbors)

Given an array where res[i] <- the neighbors of node i, return the
corresponding lightgraph.
"""
function _gen_graph(neighbors)
    res = LightGraphs.SimpleGraph(length(neighbors))

    for i in 1:length(neighbors)
        for j in 1:length(neighbors[i])
            LightGraphs.add_edge!(res, i => neighbors[i][j])
        end
    end

    return res
end

"""
    plot_faces(faces, points; plot_name, col, new_plot, individually)

Given a set of vectors corresponding to faces of a planar graph, plots the faces in the unit square.
"""
function plot_faces(faces::Set{Vector{T}}, points::Vector{Any}; plot_name::String="Faces", col::String="green3", new_plot::Bool=true, individually::Bool=false) where {T}
    if new_plot && !individually
        Plots.plot()
    end

    for (j,face) in enumerate(faces)
        v = Polyhedra.convexhull(map(i -> collect(points[i]), face)...)
        x_locations = map(i -> points[i].first, face)
        y_locations = map(i -> points[i].second, face)

        # avg_x = Statistics.mean(x_locations)
        # avg_y = Statistics.mean(y_locations)

        polygon = Polyhedra.polyhedron(
            v,
            Polyhedra.DefaultLibrary{Rational{Int64}}(GLPK.Optimizer)
        )

        if individually
            # Plot each face on a separate plot with node labels
            Plots.plot(polygon, color=col, alpha=0.9)
            Plots.plot!(x_locations, y_locations, series_annotations=([Plots.text(string(x), :center, 8, "courier") for x in face]))
            display(Plots.plot!(xlims=(-0.05,1.05), ylims=(-0.05,1.05), title="Face $j: $face"))
        else
            # Accumulate the faces on the same plot
            Plots.plot!(polygon, color=col, alpha=0.9)
            # Plots.plot!([avg_x], [avg_y], series_annotations=([Plots.text("$j", :center, 8, "courier")]))
            # display(Plots.plot!(polygon, title="Free Face # $j", xlims=(-0.05,1.05), ylims=(-0.05,1.05)))
        end
    end

    if !individually
        display(Plots.plot!(title=plot_name, xlims=(-0.05,1.05), ylims=(-0.05,1.05)))
    end
end

"""
    plot_edges(lg, points; plot_name, col, new_plot)

Given a LabeledGraph, plots its nodes and edges in the unit square.
"""
function plot_edges(lg::LabeledGraph{T}, points::Vector{Any}; plot_name::String="Edges", col::String="colorful", new_plot::Bool=true, vertices::Dict{T,T}=Dict{T,T}()) where {T}
    if new_plot
        Plots.plot()
    end

    # Need to map nodes to the point they refer to in "points"
    rev = ClutteredEnvPathOpt._reverse_labels(lg.labels)
    for edge in LightGraphs.edges(lg.graph)
        if col == "colorful"
            Plots.plot!([points[rev[edge.src]].first, points[rev[edge.dst]].first], [points[rev[edge.src]].second, points[rev[edge.dst]].second],linewidth=2)
            # display(Plots.plot!(title="Edge ($(rev[edge.src]), $(rev[edge.dst]))"))
        else
            Plots.plot!([points[rev[edge.src]].first, points[rev[edge.dst]].first], [points[rev[edge.src]].second, points[rev[edge.dst]].second], color=col,linewidth=2)
            # display(Plots.plot!(title="Edge ($(rev[edge.src]), $(rev[edge.dst]))"))
        end
    end

    # ClutteredEnvPathOpt.plot_intersections(obstacles, vertices=vertices)
    ClutteredEnvPathOpt.plot_points(points, vertices=vertices)

    display(Plots.plot!(title=plot_name, xlims=(-0.05,1.05), ylims=(-0.05,1.05), legend=false))
end

"""
    plot_field(field)

Plots the obstacles to the existing active plot.
"""
function plot_field(field)
    for i = 1:length(field)
        Plots.plot!(field[i], xlims = (-0.05,1.05), ylim = (-0.05, 1.05))
    end
end

"""
    plot_lines(field)

Plots the lines that make up the obstacles' halfspaces to the existing active
plot.
"""
function plot_lines(field)
    halfspaces = @pipe map(obstacle -> Polyhedra.hrep(obstacle).halfspaces, field) |> Iterators.flatten(_) |> collect(_)
    unique!(halfspaces)

    for h in halfspaces
        if abs(h.a[2]) != 0//1
            f = x -> (h.β - h.a[1] * x) // h.a[2]

            x = LinRange(0//1, 1//1, 11)
            y = map(f, x)
        else  # vertical lines
            x = fill(abs(h.β // h.a[1]), 11)
            y = LinRange(0//1, 1//1, 11)
        end

        Plots.plot!(x, y)
    end
end

"""
    plot_borders()

Plots the lines that make up the unit square's borders to the existing active
plot.
"""
function plot_borders()
    halfspaces = [
        Polyhedra.HalfSpace([-1//1, 0//1], 0//1),
        Polyhedra.HalfSpace([1//1, 0//1], 1//1),
        Polyhedra.HalfSpace([0//1, -1//1], 0//1),
        Polyhedra.HalfSpace([0//1, 1//1], 1//1)
    ]

    for h in halfspaces
        if abs(h.a[2]) != 0//1
            f = x -> (h.β - h.a[1] * x) // h.a[2]

            x = LinRange(0//1, 1//1, 11)
            y = map(f, x)
        else  # vertical lines
            x = fill(abs(h.β // h.a[1]), 11)
            y = LinRange(0//1, 1//1, 11)
        end

        Plots.plot!(x, y)
    end
end

"""
    plot_intersections(field)

Plots and labels the points where the lines that make up the obstacles'
halfspaces intersect to the existing active plot.
"""
function plot_intersections(field; vertices::Dict{T,T}=Dict{T,T}()) where {T}
    intersections, _, inside_quant = find_intersections(field)

    # Remove points inside obstacles
    for _ = 1:inside_quant
        pop!(intersections)
    end

    if !isempty(vertices)
        points_filtered = []
        for v in keys(vertices)
            push!(points_filtered, intersections[v])
        end
        x = map(point -> point.first, points_filtered)
        y = map(point -> point.second, points_filtered)
        Plots.scatter!(x,y, color="red", series_annotations=([Plots.text(string(x), :right, 8, "courier") for x in keys(vertices)]))
    else
        x = map(point -> point.first, intersections)
        y = map(point -> point.second, intersections)
        Plots.scatter!(x,y, color="red", series_annotations=([Plots.text(string(x), :right, 8, "courier") for x in 1:length(points)]))
    end

    # TODO: Delete this once certain it works
    # if !isempty(vertices)
    #     intersections = @pipe filter(v -> v in keys(vertices), 1:length(intersections)) |> intersections[_]
    # end

    # x = map(point -> point[1], intersections)
    # y = map(point -> point[2], intersections)

    # # Plots.scatter!(x,y, color="red3", series_annotations=([Plots.text(string(x), :right, 8, "courier") for x in 1:length(x)]))
    # Plots.scatter!(x,y, color="red", series_annotations=([Plots.text(string(x), :right, 8, "courier") for x in 1:length(x)]))
end

"""
    plot_points(points; vertices)

Plots and labels points given.
"""
function plot_points(points::Vector{Any}; vertices::Dict{Int,Int}=Dict{Int,Int}())

    Plots.plot!(legend=false)

    if !isempty(vertices)
        points_filtered = []
        for v in keys(vertices)
            push!(points_filtered, points[v])
        end
        x = map(point -> point.first, points_filtered)
        y = map(point -> point.second, points_filtered)
        Plots.scatter!(x,y, color="red", series_annotations=([Plots.text(string(x), :right, 8, "courier") for x in keys(vertices)]))
    else
        x = map(point -> point.first, points)
        y = map(point -> point.second, points)
        Plots.scatter!(x,y, color="red", series_annotations=([Plots.text(string(x), :right, 8, "courier") for x in 1:length(points)]))
    end

end

"""
    plot_new(n, name; custom, seed, save_image)

Generates a new obstacle course with n obstacles. If custom = true, obstacles will be
created from an array containing the points for each obstacle. If save_image = true,
the course will be plotted and saved to an image.
"""
function plot_new(n::Int, name::String; custom::Bool=false, seed::Int=1, save_image::Bool=false, partition::String="CDT")
    if custom
        # obs = gen_field(num_obstacles::Int; custom::Bool=false, points_it::Set{Vector{T}}=Set{Vector{Rational}}(), seed::Int=1) where {T}
    else
        obs = gen_field(n, seed = seed)
    end

    if save_image
        Plots.plot()
        plot_field(obs)
        points = ClutteredEnvPathOpt.find_points(obs)
        ClutteredEnvPathOpt.plot_points(points)
        #plot_lines(obs)
        #plot_borders()
        #plot_intersections(obs)
        Plots.png(name)
        display(Plots.plot!(title="Field"))
    end

    if partition == "CDT"
        return construct_graph_delaunay(obs)
    else
        return construct_graph(obs)
    end
end


text:
# custom commands (lx_* functions; create a Utils module and put stuff in there)

# custom environments (env_* functions; )

ext: jl


text:
using SparseQuadratureGrids
using Base.Test

# write your own tests here

Base.var(::Type{GenzKeister}) = 1/2
Base.var(::Type{KronrodPatterson}) = 1/3

p = 40, n = 1000
Am = randn(p, n);
Bsv = [StaticArrays.SVector{p}(randn(p)) for i ∈ 1:n];
Ab = SparseQuadratureGrids.MatrixVecSVec(Am, Val{p});
Bb = SparseQuadratureGrids.MatrixVecSVec(Bsv);

function testSplitWeights(g::SplitWeights{p,q}, ec::Int) where {p,q}
  @testset begin
    @test sum(g.weight_sum) ≈ 1.0
    @test g.total_count == ec
    @test minimum(g.nodes_positive) < minimum(g.nodes_negative)
    @test maximum(g.nodes_positive) > maximum(g.nodes_negative)
    @test all(isapprox.(g.nodes_positive * g.weights_positive, 0.0, atol = 1e-14))
    @test all(isapprox.(g.nodes_negative * g.weights_negative, 0.0, atol = 1e-14))
    @test all(isapprox.(g.nodes_positive .^ 2 * g.weights_positive .+ g.nodes_negative .^ 2 * g.weights_negative, var(q)))
  end
end

dsgk = [1, 3, 3, 9, 9, 9, 9, 19, 19, 19, 19, 19, 19, 19, 19, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 35, 103]
dskp = [1, 3, 3, 7, 7, 7, 15, 15, 15, 15, 15, 15, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 63]
testSplitWeights(SplitWeights(10,6,GenzKeister), 115645)
testSplitWeights(SplitWeights(4,length(dsgk),GenzKeister,dsgk), 253649)

testSplitWeights(SplitWeights(8,6,KronrodPatterson), 31745)
testSplitWeights(SplitWeights(5,length(dskp),KronrodPatterson,dskp), 1297727)


text:
using Test
using Libxc

@testset "Version" begin
    @test isa(Libxc.xc_version(), VersionNumber)
    @test isa(Libxc.xc_version_string(), String)
end

@testset "FuncReferType" begin
    refer = Libxc.FuncReferType(pointer("ref"), pointer("doi"), pointer("bibtex"))
    p = Ref{Libxc.FuncReferType}(refer)
    @test unsafe_string(Libxc.xc_func_reference_get_ref(p)) == "ref"
    @test unsafe_string(Libxc.xc_func_reference_get_doi(p)) == "doi"
    @test unsafe_string(Libxc.xc_func_reference_get_bibtex(p)) == "bibtex"
end

@testset "XCFuncInfoType" begin
    reff = Libxc.xc_func_alloc()
    Libxc.xc_func_init(reff, 1, 1)
    info_ref = Libxc.xc_func_get_info(reff)
    @test Libxc.xc_func_info_get_number(info_ref) == 1
    @test Libxc.xc_func_info_get_kind(info_ref) == 0
    @test Libxc.xc_func_info_get_name(info_ref) == "Slater exchange"
    @test Libxc.xc_func_info_get_family(info_ref) == 1
    @test Libxc.xc_func_info_get_flags(info_ref) == 143
    Libxc.xc_func_end(reff)
    Libxc.xc_func_free(reff)
end

@testset "XCFuncType xc calc" begin
    rho = [0.1, 0.2, 0.3, 0.4, 0.5]
    sigma = [0.2, 0.3, 0.4, 0.5, 0.6]
    result = zeros(Float64, 5)

    ptr = Libxc.xc_func_alloc()
    Libxc.xc_func_init(ptr, Libxc.LDA_X, 1)
    Libxc.xc_lda_exc!(ptr, 5, rho, result)
    @test result ≈ [-0.342809, -0.431912, -0.494416, -0.544175, -0.586194] atol=1e-5
    Libxc.xc_func_end(ptr)
    Libxc.xc_func_free(ptr)

    ptr = Libxc.xc_func_alloc()
    Libxc.xc_func_init(ptr, Libxc.GGA_X_PBE, 1)
    Libxc.xc_gga_exc!(ptr, 5, rho, sigma, result)
    @test result ≈ [-0.452598, -0.478878, -0.520674, -0.561428, -0.598661] atol=1e-5
    Libxc.xc_func_end(ptr)
    Libxc.xc_func_free(ptr)
end

@testset "High-level interface" begin
    rho = [0.1, 0.2, 0.3, 0.4, 0.5]
    sigma = [0.2, 0.3, 0.4, 0.5, 0.6]
    result = zeros(Float64, 5)

    # LDA
    func = Libxc.Functional(:lda_x)
    @test func.identifier == :lda_x
    @test func.number == 1
    @test func.family == Libxc.family_lda
    @test func.n_spin == 1

    Libxc.evaluate_lda!(func, rho, E=result)
    @test result ≈ [-0.342809, -0.431912, -0.494416, -0.544175, -0.586194] atol=1e-5

    # GGA
    func = Libxc.Functional(:gga_x_pbe)
    @test func.identifier == :gga_x_pbe
    @test func.number == 101
    @test func.family == Libxc.family_gga
    @test func.n_spin == 1

    Libxc.evaluate_gga!(func, rho, sigma, E=result)
    @test result ≈ [-0.452598, -0.478878, -0.520674, -0.561428, -0.598661] atol=1e-5
end


text:
#=
AirFoilCoordinatesArraysActivity.jl
Jacob Child
April 23rd, 2022

Equations Used
Zu(Zupper) = zbar + t/2
Zl(Zlower) = zbar -t/2
t = 10m(.2969sqrt(x) - .1620x - .3537x^2 + .2843x^3 - .1015x^4)
if X<=p
        zbar = c*(2*p*x-x^2)/(p^2)
    elseif x > p
        zbar = c*(1-2*p+2*p*x-x^2)/((1-p)^2)
    end

Pseudocode- Take in NACA 4- Series Airfoil numbers- ie max camber (c),
    max camber position (p), thickness (t), and an array (x) to plot at.
    For example NACA 2412 = c: .02, p: .04, t: .12 (all in meters I believe)
    Do that by asking the user for c,p,t, then creating x (hardcode it for simplicity)
    Use the functions previously written to calculate zu and zl.

    Goal = clean code, lots of well commented functions
=#

"""
user_input Function to take in the 4 NACA numbers
"""
function user_input()
    println("Input the 4 NACA numbers ie 2412 \nNACA ")
    input1 = readline()
    println("\nDo you want Cosine Spacing? 1 = Yes, 0 = No\n")
    input2 = readline()
    if input2 == "1"
        println("\nHow many points do you want?\n")
        stepsize = pi / parse(Int64,readline() )
        #println("\n", stepsize)
        theta = collect(0.0:stepsize:pi)    #creates numbers from 0 to pi
        x = (1/2) * (1 .- cos.(theta))  #should give proper Spacing, assuming c (chord length) = 1

    else
        println("\nNot using Cosine Spacing")
        x = collect(0:0.1:1)   #10 points

    end

    c = parse(Int64, input1[1])
    c = c/100  #converts the chord to the proper units   #can do c//100 and that leaves it in fraction form
    p = parse(Int64, input1[2])
    p = p/10    #converts the max camber position to the proper units
    t = parse(Int64, input1[3:4])
    t = t/100
    #println("\n c:", c, "\n p:", p, "\n t:", t)    #prints out the actual NACA values

    return c, p, t, input1, x
end

"""
thickness_calculations(x,p)
Takes in the x coordinates, and max thickness (p), calculates thickness, and ouputs the array

"""
function thickness_calculations(x,t)
    m = t
    TArray = 10*m*(.2969*sqrt(x) - .1260*x - .3537*x^2 + .2843*x^3 - .1015*x^4)
    return TArray
end

"""
zbar_calculations(c,p,x)
Takes in c, p, and x values and outputs zbar
"""
function zbar_calculations(c,p,x)
    zbar = Float64[] #makes zbar an array of zeros the length of x, hopefully?
    for i in x
        if i <= p
             zbarInput = c*(2*p*i - i^2)/(p^2)
             push!(zbar,zbarInput)

        elseif i > p
             zbarInput = c*(1 - 2*p + 2*p*i - i^2)/((1-p)^2)
             push!(zbar,zbarInput)

        else
             zbar = "you wenti"
             break

        end

    end
    return zbar
end

"""
zupper_lower(zbar, TArray)
calculates zu and zl and returns both
"""
function zupper_lower(zbar, TArray)

    zu = zbar .+ TArray/2
    zl = zbar .- TArray/2

    return zu, zl
end

#Actual code- See how it goes

c, p, t, fromuser, x = user_input()  #function to get user input
#why is collect better than x = range(0:0.010:1)?
TArray = thickness_calculations.(x, t)   #calculates the thicknesses
zbarn = zbar_calculations(c, p, x) #calculate and returns zbar
zu, zl = zupper_lower(zbarn, TArray) #calculates and returns zupper and z lower

#outputs

#println("\n \n", zu, "\n \n", zl) #for vs code

using DataFrames #only seems to work in the Julia REPL

DataFrame()

Data = DataFrame("X" => x, "Zupper"=> zu, "Zlower" => zl, "Thickness" => TArray, "Zbar" => zbarn)

display(Data)

#look into aspect ratio for plots
 using Plots
 plt = plot(x,zu)   #saving the plot to an element
 plot!(x,zl)    #adding/updating to that plot
 #Plot settings
 plot!(xlims = (-.1, 1.1), title = ("NACA $fromuser"), xlabel = "Position", ylabel = "Position")
 plot!(aspectratio = :equal)

 display(plt)   #displaying the plot by calling the saved element


text:
"""
    PyIterable{T=Py}(x)

This object iterates over iterable Python object `x`, yielding values of type `T`.
"""
struct PyIterable{T}
    py :: Py
    PyIterable{T}(::Val{:new}, py::Py) where {T} = new{T}(py)
end
export PyIterable

PyIterable{T}(x) where {T} = PyIterable{T}(Val(:new), Py(x))
PyIterable(x) = PyIterable{Py}(x)

ispy(x::PyIterable) = true
getpy(x::PyIterable) = x.py
pydel!(x::PyIterable) = pydel!(x.py)

Base.IteratorSize(::Type{PyIterable{T}}) where {T} = Base.SizeUnknown()
Base.eltype(::Type{PyIterable{T}}) where {T} = T

function Base.iterate(x::PyIterable{T}, it::Py=pyiter(x)) where {T}
    y = pynext(it)
    if ispynull(y)
        pydel!(it)
        return nothing
    else
        return (pyconvert_and_del(T, y), it)
    end
end

pyconvert_rule_iterable(::Type{T}, x::Py, ::Type{PyIterable{V}}=Utils._type_ub(T)) where {T<:PyIterable,V} =
    if PyIterable{Py} <: T
        pyconvert_return(PyIterable{Py}(x))
    else
        pyconvert_return(PyIterable{V}(x))
    end


text:
abstract type TullyModel <: DiabaticModel end

NQCModels.ndofs(::TullyModel) = 1
NQCModels.nstates(::TullyModel) = 2

"""
    TullyModelOne(a=0.01, b=1.6, c=0.005, d=1.0)

Tully's simple avoided crossing model from [J. Chem. Phys. 93, 1061 (1990)](https://doi.org/10.1063/1.459170).
"""
Parameters.@with_kw struct TullyModelOne{A,B,C,D} <: TullyModel
    a::A = 0.01
    b::B = 1.6
    c::C = 0.005
    d::D = 1.0
end

function NQCModels.potential(model::TullyModelOne, q::Real)
    Parameters.@unpack a, b, c, d = model
    if q > 0
        V11 = a * (1 - exp(-b*q))
    else
        V11 = -a * (1 - exp(b*q))
    end
    V22 = -V11
    V12 = c * exp(-d*q^2) # sets both off-diagonals as V is Hermitian
    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative(model::TullyModelOne, q::Real)
    Parameters.@unpack a, b, c, d = model
    D11 = a * b * exp(-b * abs(q))
    D22 = -D11
    D12 = -2 * c * d * q * exp(-d*q^2)
    return Hermitian(SMatrix{2,2}(D11, D12, D12, D22))
end

"""
    TullyModelTwo(a=0.1, b=0.28, c=0.015, d=0.06, e=0.05)

Tully's dual avoided crossing model from [J. Chem. Phys. 93, 1061 (1990)](https://doi.org/10.1063/1.459170).
"""
Parameters.@with_kw struct TullyModelTwo{A,B,C,D,E} <: TullyModel
    a::A = 0.1
    b::B = 0.28
    c::C = 0.015
    d::D = 0.06
    e::E = 0.05
end

function NQCModels.potential(model::TullyModelTwo, q::Real)
    Parameters.@unpack a, b, c, d, e = model
    V11 = 0
    V22 = -a*exp(-b*q^2) + e
    V12 = c * exp(-d*q^2)
    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative(model::TullyModelTwo, q::Real)
    Parameters.@unpack a, b, c, d  = model
    D11 = 0
    D22 = 2*a*b*q*exp(-b*q^2)
    D12 = -2*c*d*q*exp(-d*q^2)

    Hermitian(SMatrix{2,2}(D11, D12, D12, D22))
end

"""
    TullyModelThree(a=0.0006, b=0.1, c=0.9)

Tully's extended coupling with reflection model from [J. Chem. Phys. 93, 1061 (1990)](https://doi.org/10.1063/1.459170).
"""
Parameters.@with_kw struct TullyModelThree{A,B,C} <: TullyModel
    a::A = 0.0006
    b::B = 0.1
    c::C = 0.9
end

function NQCModels.potential(model::TullyModelThree, q::Real)
    Parameters.@unpack a, b, c = model
    V11 = a
    V22 = -a
    if q > 0
        V12 = b * (2 - exp(-c*q))
    else
        V12 = b * exp(c*q)
    end
    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative(model::TullyModelThree, R::Real)
    Parameters.@unpack a, b, c = model
    q = R[1]
    D11 = 0
    D22 = 0
    if q > 0
        D12 = b * c * exp(-c*q)
    else
        D12 = b * c * exp(c*q)
    end
    return Hermitian(SMatrix{2,2}(D11, D12, D12, D22))
end


text:
"""
    assemble_flux_matrices(points, facets)

Assemble flux matrix (`Q`) for each compartment and boundary.
"""
function assemble_flux_matrices(points, facets)

    # Sizes
    ncompartment, nboundary = size(facets)

    flux_matrices = [
        assemble_flux_matrix(facets[icmpt, iboundary]', points[icmpt]') for
        icmpt = 1:ncompartment, iboundary = 1:nboundary
    ]

    flux_matrices
end


text:
module REALMSupport

using Images, ImageSegmentation, ImageTransformations, Optim, ImageMorphology
using Statistics, OffsetArrays
using FileIO, Printf, ProgressMeter, RegisterMismatch, Interpolations
using RegisterMismatch.RegisterCore

export grid_resolution, gauss_line_fit, img_save, beads_segment
export centerofmass, findmid, translate, translate_optim
include("grid_fun.jl")
include("gauss_fit.jl")
include("img_save.jl")
include("COM.jl")
include("translation.jl")
include("segmentbeads.jl")

end


text:
## To use this test routine, execute it
## from the REPL via
##    include(joinpath(Pkg.dir("microcoverage"), "test_microcoverage.jl"))
## It should display some messages on the
## console, and then a file called
## test_microcoverage_ex.jl.mcov should
## be generated that should be the
## same as
##   test_microcoverage_ex.jl.mcov.correct

include(joinpath(Pkg.dir("microcoverage"), "microcoverage.jl"))
using microcoverage
begintrack(joinpath(Pkg.dir("microcoverage"), "test_microcoverage_ex.jl"))
include(joinpath(Pkg.dir("microcoverage"), "test_microcoverage_ex.jl"))
test_microcoverage_ex.runfuncs()
endtrack(joinpath(Pkg.dir("microcoverage"), "test_microcoverage_ex.jl"))

max_stars_repo_name: StephenVavasis/microcoverage

max_issues_repo_name: StephenVavasis/microcoverage

max_forks_repo_name: StephenVavasis/microcoverage


text:
module BallisticLimit
using ...Types
export v_50
function equivalent_steel_cylinder(mat_p::MaterialPenetrator{WF}, mat_t::MetallicParts{ST})
    if mat_t.material.bhn<=380
        return true
    else
        return false
    end

end

function equivalent_steel_cylinder(mat_p::MaterialPenetrator{WF}, mat_t::MetallicParts{AL})
    return true
end

function equivalent_steel_cylinder(mat_p::MaterialPenetrator{WF}, mat_t::MetallicParts{TI})
    return true
end

function equivalent_steel_cylinder(mat_p::MaterialPenetrator{TA}, mat_t::MetallicParts{ST})
    return true
end

function equivalent_steel_cylinder(mat_p::MaterialPenetrator{DU}, mat_t::MetallicParts{ST})
    return true
end

function equivalent_steel_cylinder(mat_p::MaterialPenetrator{DU}, mat_t::MetallicParts{TI})
    return true
end

function normalised(penetration::Impact, proj::Projectile2D, tar_mat::AbstractMaterial)

    Ap_e = effective_area(proj,tar_mat)
    L,D,Ap_f, Ap_nf = equivalent_cylinder(proj.type, proj)
    X = penetration.los/(4*Ap_e/pi)^(1/2)
    Y = ((X*D/penetration.los)^3*(L/D)*(proj.mat.prop.density/ρ_st))^(1/2)
    flag = false
    #flag = equivalent_steel_cylinder(proj.mat, tar_mat)
    if isa(proj.mat,MaterialPenetrator{WF}) && isa(tar_mat,MetallicParts{ST})
        if proj.mat.type.bhn <= 380
            flag = true
        end
    end

    if isa(proj.mat,MaterialPenetrator{WF}) && isa(tar_mat,MetallicParts{AL})
        flag = true
    end

    if isa(proj.mat,MaterialPenetrator{WF}) && isa(tar_mat,MetallicParts{TI})
        flag = true
    end

    if isa(proj.mat,MaterialPenetrator{TA}) && isa(tar_mat,MetallicParts{ST})
        flag = true
    end

    if isa(proj.mat,MaterialPenetrator{DU}) && isa(tar_mat,MetallicParts{ST})
        flag = true
    end

    if isa(proj.mat,MaterialPenetrator{DU}) && isa(tar_mat,MetallicParts{TI})
        flag = true
    end

    if flag
        Z =Y
    else
        Z = Y*(ρ_st/proj.mat.prop.density)^(1/2)
    end
    return X, Y, Z
end

function coefficient(tar_mat::MetallicParts{ST}, proj::Projectile2D)
    if isa(proj.mat,MaterialPenetrator{ST})||isa(proj.mat,MaterialPenetrator{AL})
        C = 4254
        K = -538
        J = 5064
        b = 0.61
    elseif isa(proj.mat,MaterialPenetrator{WF})||isa(proj.mat,MaterialPenetrator{DU})
        C = 4599
        K = -594
        J = 5343
        b = 0.61
    elseif isa(proj.mat,MaterialPenetrator{TA})
        C = 2695
        K = 136
        J = 2916
        b = 1.24
    end
    return C, K, J, b
end

function coefficient(tar_mat::MetallicParts{AL}, proj::Projectile2D)
    if isa(proj.mat,MaterialPenetrator{WF})
        C = 946
        K = 0
        J = 1300
        b = 0.86
    else
        C = 3.87 * tar_mat.type.bhn + 295
        K = 2.2 * tar_mat.typ.bhn + 715
        J = -21.6 * tar_mat.type.bhn + 7282
        b = 1.75
    end
    return C, K, J, b

end

function coefficient(tar_mat::MetallicParts{TI}, proj::Projectile2D)
    if isa(proj.mat,MaterialPenetrator{WF})||isa(proj.mat,MaterialPenetrator{DU})
        C = 1495
        K = 200
        J = 3490
        b = 1.0
    else
        C = 2440
        K = 410
        J = 6275
        b = 1.05
    end
    return C, K, J, b

end

function coefficient_AP(tar_mat::MetallicParts{AL}, proj::Projectile2D)
    C_s = 1795
    bs = 0.61
    ps=1.0
    return C_s, bs, ps

end

function coefficient_AP(tar_mat::MetallicParts{TI}, proj::Projectile2D)
    C_s = 1437
    bs = 0.75
    ps=1.0
    return C_s, bs, ps

end

function coefficient_AP(tar_mat::MetallicParts{ST}, proj::Projectile2D)
    if isa(proj.mat,MaterialPenetrator{DU}) || isa(proj.mat,MaterialPenetrator{WF}) || isa(proj.mat,MaterialPenetrator{TA})
        C_s = 2763
        bs = 0.66
        ps = 1.0
    else
        C_s = 3143
        bs = 0.60
        ps = 1.0
    end

    return C_s, bs, ps

end
function Vfactor(proj::Projectile2D, tar_mat::AbstractMaterial, penetration::Impact)
    Ap_e = effective_area(proj,tar_mat)
    γ_θ = tar_mat.prop.density*penetration.los*Ap_e/(proj.mass*cosd(penetration.obl))
    γ_θr =  tar_mat.prop.density*penetration.los*Ap_e/(proj.mass*cosd(θ_ref))
    V_factor = (cosd(θ_ref)/cosd(penetration.obl))*sqrt((exp(γ_θ)-1)/(exp(γ_θr)-1))
    return V_factor
end

function ballistic_limit(penetration::Impact, proj::Projectile2D, tar_mat::MetallicParts)

    θ_l = penetration.obl
    if penetration.obl > θ_ref
        θ_l = θ_ref
    end
    X, Y, Z = normalised(penetration, proj, tar_mat)
    C, K, J, b = coefficient(tar_mat, proj)
    if X>0.1
        V_50n = (C*X^b+K)/Y
    else
        V_50n = J*X/Y
    end
    if isa(tar_mat,MetallicParts{ST})
        C_st = -6.37*10^-6*tar_mat.prop.bhn^2+0.00373*tar_mat.prop.bhn+0.456
        #println(V_50n)
        V_50n = C_st*V_50n
        #println(V_50n)
    end
    #if typeof(tar_mat)==GE
    #    V_50 = (0.985*θ_l^2-0.74*θ_l+1.0)*V_50n
    #elseif X<=2
    #    V_50 = (1-((0.197-0.0486*X^2.45)*θ_l^2+(0.0745*X^2-0.135*X)*θ_l))*(V_50n/cosd(θ_l))
    #else
        #println(V_50n)
        V_50 = V_50n/cosd(θ_l)
    #end
    if penetration.obl > θ_ref
        V_factor = Vfactor(proj, tar_mat, penetration)
        #println(V_50)
        V_50 = V_50*V_factor
    end
    #println(V_50)
    return V_50
end

function nose_type(AP::AP)
    if (AP.nose_diameter/AP.base_diameter<=0.5)&&(AP.nose_length/(AP.base_diameter-AP.nose_diameter)>=1.0)
        NSHP = false
    else
        NSHP = true
    end
    return NSHP
end

function attack_type(proj::Projectile2D, tar_mat::AbstractMaterial)
    NSHP = nose_type(proj.type)
    if NSHP || (proj.yaw>10)||(isa(tar_mat,MetallicParts{ST}) && proj.mat.prop.bhn<560)
        ATK = false #blunt attack
    else
        ATK = true #sharp attack
    end
    return ATK
end

function ballistic_limit_sharp(penetration::Impact, proj::Projectile2D, tar_mat::MetallicParts, proj_type::AbstractProjectile)
    X, Y, Z = normalised(penetration, proj, tar_mat)
    L,D,Ap_f, Ap_nf = equivalent_cylinder(proj.type, proj)
    C_s, bs, ps = coefficient_AP(tar_mat, proj)
    V_50n = C_s*((ρ_st/proj.mat.prop.density)*(X/(L/D)))^bs*(1-exp(-3.2*sqrt(X)))
    θ = penetration.obl
    if isa(tar_mat,MetallicParts{GE})
        V_50 = V_50n*(1.640*(deg2rad(θ))^2-0.72*deg2rad(θ)+1.0)
    else
        V_50 = V_50n/(cosd(θ))^ps
    end
    if θ > θ_ref
        V_factor = Vfactor(proj, tar_mat, penetration)
        V_50 = V_50 * V_factor
    end

    return V_50

end

function ballistic_limit_sharp(penetration::Impact, proj::Projectile2D, tar_mat::MaterialUnkown, proj_type::AbstractProjectile)
    return 0
end
function ballistic_limit(penetration::Impact, proj::Projectile2D, tar_mat::MaterialUnkown, proj_type::AbstractProjectile)
    return 0
end

#function ballistic_limit_sharp(penetration::Impact, proj::Projectile, tar_mat::Crew, proj_type::AbstractProjectile)
#    return 0
#end
#function ballistic_limit(penetration::Impact, proj::Projectile, tar_mat::Crew, proj_type::AbstractProjectile)
#    return 0
#end

function v_50(proj::AbstractPenetrator, tar_mat::AbstractMaterial, penetration::Impact)
    if typeof(proj.type) <: AbstractFragment
        v50 = ballistic_limit(penetration, proj, tar_mat)
        #println(v50)
    elseif typeof(proj.type) <: AbstractProjectile
        sharp = attack_type(proj, tar_mat)
        if sharp
            v50 = ballistic_limit_sharp(penetration, proj, tar_mat, proj.type)
        else
            v50 = ballistic_limit(penetration, proj, tar_mat)
        end
    end
    #println(v50)
    v50 = v50*0.3048  #in m/s
    return v50
end
end


text:
function __init__()
    global DEFAULT_RESOURCE = Ref{AbstractResource}(CPU1())

    TRAIT_FUNCTION_GIVEN_NAME[:measure]      = is_measure
    TRAIT_FUNCTION_GIVEN_NAME[:measure_type] = is_measure_type

    MLJModelInterface.set_interface_mode(MLJModelInterface.FullInterface())
end


text:
### Latent models ####
include("prior.jl")
include("posterior.jl")

## Exact Gaussian Process
struct LatentGP{T,Tpr<:GPPrior,Tpo<:Posterior{T},O} <: AbstractLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
    opt::O
end

function LatentGP(T::DataType, dim::Int, kernel::Kernel, mean::PriorMean, opt)
    return LatentGP(
        GPPrior(deepcopy(kernel), deepcopy(mean)),
        Posterior(dim, zeros(T, dim), cholesky(Matrix{T}(I(dim)))),
        deepcopy(opt),
    )
end

@traitimpl IsFull{LatentGP}

## AbstractVarLatent

abstract type AbstractVarLatent{T,Tpr,Tpo} <: AbstractLatent{T,Tpr,Tpo} end

## Variational Gaussian Process
mutable struct VarLatent{T,Tpr<:GPPrior,Tpo<:VarPosterior{T},O} <:
               AbstractVarLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
    opt::O
end

function VarLatent(T::DataType, dim::Int, kernel::Kernel, mean::PriorMean, opt)
    return VarLatent(
        GPPrior(deepcopy(kernel), deepcopy(mean)), VarPosterior{T}(dim), deepcopy(opt)
    )
end

@traitimpl IsFull{VarLatent}

## Sparse Variational Gaussian Process

mutable struct SparseVarLatent{
    T,Tpr<:GPPrior,Tpo<:VarPosterior{T},Topt,TZ<:AbstractVector,TZopt
} <: AbstractVarLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
    Z::TZ
    opt::Topt
    Zopt::TZopt
end

function SparseVarLatent(
    T::DataType,
    Z::AbstractVector,
    kernel::Kernel,
    mean::PriorMean,
    opt=nothing,
    Zopt=nothing,
)
    dim = length(Z)
    return SparseVarLatent(
        GPPrior(deepcopy(kernel), deepcopy(mean)),
        VarPosterior{T}(dim),
        deepcopy(Z),
        deepcopy(opt),
        deepcopy(Zopt),
    )
end

@traitimpl IsSparse{SparseVarLatent}

## Monte-Carlo Gaussian Process

struct SampledLatent{T,Tpr<:GPPrior,Tpo<:SampledPosterior{T}} <: AbstractLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
end

function SampledLatent(T::DataType, dim::Int, kernel::Kernel, mean::PriorMean)
    return SampledLatent(
        GPPrior(deepcopy(kernel), deepcopy(mean)),
        SampledPosterior(dim, zeros(T, dim), Symmetric(Matrix{T}(I(dim)))),
    )
end

@traitimpl IsFull{SampledLatent}

## Online Sparse Variational Process

mutable struct OnlineVarLatent{
    T,
    Tpr<:GPPrior,
    Tpo<:AbstractVarPosterior{T},
    Topt,
    TZ<:AbstractVector,
    TZalg<:InducingPoints.OnIPSA,
    TZopt,
} <: AbstractVarLatent{T,Tpo,Tpr}
    prior::Tpr
    post::Tpo
    Z::TZ
    Zₐ::TZ
    Zalg::TZalg
    Zupdated::Bool
    opt::Topt
    Zopt::TZopt
end

function OnlineVarLatent(
    T::DataType,
    dim::Int,
    Z::AbstractVector,
    Zalg::InducingPoints.OnIPSA,
    kernel::Kernel,
    mean::PriorMean,
    opt=nothing,
    Zopt=nothing,
)
    return OnlineVarLatent(
        GPPrior(deepcopy(kernel), deepcopy(mean)),
        OnlineVarPosterior{T}(dim),
        Z,
        deepcopy(Z),
        Zalg,
        false,
        deepcopy(opt),
        deepcopy(Zopt),
    )
end

@traitimpl IsSparse{OnlineVarLatent}

## Variational Student-T Process

mutable struct TVarLatent{T<:Real,Tpr<:TPrior,Tpo<:VarPosterior{T},O} <:
               AbstractVarLatent{T,Tpr,Tpo}
    prior::Tpr
    post::Tpo
    opt::O
end

function TVarLatent(T::DataType, ν::Real, dim::Int, kernel::Kernel, mean::PriorMean, opt)
    return TVarLatent(
        TPrior(deepcopy(kernel), deepcopy(mean), ν, rand(T), rand(T)),
        VarPosterior{T}(dim),
        deepcopy(opt),
    )
end

@traitimpl IsFull{TVarLatent}

### Functions

prior(gp::AbstractLatent) = gp.prior
kernel(gp::AbstractLatent) = kernel(prior(gp))
setkernel!(gp::AbstractLatent, kernel::Kernel) = setkernel!(prior(gp), kernel)
pr_mean(gp::AbstractLatent) = mean(prior(gp))
pr_mean(gp::AbstractLatent, X::AbstractVector) = mean(prior(gp), X)
setpr_mean!(gp::AbstractLatent, μ₀::PriorMean) = setmean!(prior(gp), μ₀)

posterior(gp::AbstractLatent) = gp.post
Distributions.dim(gp::AbstractLatent) = dim(posterior(gp))
Distributions.mean(gp::AbstractLatent) = mean(posterior(gp))
Distributions.cov(gp::AbstractLatent) = cov(posterior(gp))
Distributions.var(gp::AbstractLatent) = var(posterior(gp))
nat1(gp::AbstractVarLatent) = nat1(posterior(gp))
nat2(gp::AbstractVarLatent) = nat2(posterior(gp))

@traitfn function mean_f(gp::T, kernel_matrices) where {T <: AbstractLatent; IsFull{T}}
    return mean_f(mean(gp))
end
@traitfn function mean_f(gp::T, kernel_matrices) where {T <: AbstractLatent; !IsFull{T}}
    return mean_f(mean(gp), kernel_matrices.κ)
end

mean_f(μ::AbstractVector) = μ
mean_f(μ::AbstractVector, κ::AbstractMatrix) = κ * μ

@traitfn function var_f(gp::T, kernel_matrices) where {T <: AbstractLatent; IsFull{T}}
    return var_f(cov(gp))
end
@traitfn function var_f(gp::T, kernel_matrices) where {T <: AbstractLatent; !IsFull{T}}
    return var_f(cov(gp), kernel_matrices.κ, kernel_matrices.K̃)
end

var_f(Σ::AbstractMatrix) = diag(Σ)
var_f(Σ::AbstractMatrix, κ::AbstractMatrix, K̃::AbstractVector) = diag_ABt(κ * Σ, κ) + K̃

Zview(gp::SparseVarLatent) = gp.Z
Zview(gp::OnlineVarLatent) = gp.Z

setZ!(gp::AbstractLatent, Z::AbstractVector) = gp.Z = Z

opt(gp::AbstractLatent) = gp.opt
Zopt(::AbstractLatent) = nothing
Zopt(gp::SparseVarLatent) = gp.Zopt
Zopt(gp::OnlineVarLatent) = gp.Zopt

function compute_K(gp::AbstractLatent, X::AbstractVector, jitt::Real)
    return cholesky(kernelmatrix(kernel(gp), X) + jitt * I)
end

function compute_K(gp::AbstractLatent, jitt::Real)
    return cholesky(kernelmatrix(kernel(gp), gp.Z) + jitt * I)
end

function compute_κ(gp::SparseVarLatent, X::AbstractVector, K, jitt::Real)
    Knm = kernelmatrix(kernel(gp), X, gp.Z)
    κ = copy(Knm / K)
    K̃ = kernelmatrix_diag(kernel(gp), X) .+ jitt - diag_ABt(κ, Knm)
    all(K̃ .> 0) || error("K̃ has negative values")
    return (; Knm, κ, K̃)
end

function compute_κ(gp::OnlineVarLatent{T}, X::AbstractVector, K, jitt::Real) where {T}
    k = dim(gp)
    # Covariance with the model at t-1
    if isempty(gp.Zₐ)
        Kab = zeros(T, k, k)
        κₐ = I
        K̃ₐ = zero(Kab)
    else
        Kab = kernelmatrix(kernel(gp), gp.Zₐ, gp.Z)
        κₐ = Kab / K
        Kₐ = Symmetric(kernelmatrix(kernel(gp), gp.Zₐ) + jitt * I)
        K̃ₐ = Kₐ - κₐ * transpose(Kab)
    end

    # Covariance with a new batch
    Knm = kernelmatrix(kernel(gp), X, gp.Z)
    κ = Knm / K
    K̃ = kernelmatrix_diag(kernel(gp), X) .+ jitt - diag_ABt(κ, Knm)
    all(K̃ .> 0) || error("K̃ has negative values")
    return (; Kab, κₐ, K̃ₐ, Knm, κ, K̃)
end


text:
using JLD2
using NCDatasets
using OrderedCollections
using Test
using CLIMA.Writers

@testset "Writers" begin
    odims = OrderedDict(
        "x" => collect(1:5),
        "y" => collect(1:5),
        "z" => collect(1010:10:1050),
    )
    ovars = OrderedDict("v1" => rand(5, 5, 5), "v2" => rand(5, 5, 5))
    jfn, _ = mktemp()
    jfull = full_name(JLD2Writer(), jfn)
    write_data(JLD2Writer(), jfn, odims, ovars)
    nfn, _ = mktemp()
    nfull = full_name(NetCDFWriter(), nfn)
    write_data(NetCDFWriter(), nfn, odims, ovars)

    jldopen(jfull, "r") do jds
        @test get(jds, "dim_1", "foo") == "x"
        @test get(jds, "dim_2", "foo") == "y"
        @test get(jds, "dim_3", "foo") == "z"
        @test get(jds, "dim_4", "foo") == "foo"
        @test jds["x"] == odims["x"]
        @test jds["y"] == odims["y"]
        @test jds["z"] == odims["z"]
        @test get(jds, "v1", []) == ovars["v1"]
        @test get(jds, "v2", []) == ovars["v2"]
    end

    NCDataset(nfull, "r") do nds
        @test nds["x"] == odims["x"]
        @test nds["y"] == odims["y"]
        @test nds["z"] == odims["z"]
        @test try
            nds["a"] == ones(5)
            false
        catch e
            true
        end
        @test nds["v1"] == ovars["v1"]
        @test nds["v2"] == ovars["v2"]
    end
end

max_stars_repo_name: trontrytel/CLIMA

max_issues_repo_name: trontrytel/CLIMA

max_forks_repo_name: trontrytel/CLIMA


text:
################################################################################
#                              ACCESSER METHODS
################################################################################
"""
    collection_sets(set::AbstractInfiniteSet)

Return the array of sets associated with a `CollectionSet`. Error if the input
set is not a `CollectionSet`.
"""
function collection_sets(set::AbstractInfiniteSet)
    type = typeof(set)
    error("`collection_sets` not defined for infinite set of type $type.")
end
function collection_sets(set::CollectionSet{S}
                         )::Vector{S} where {S <: InfiniteScalarSet}
    return set.sets
end

################################################################################
#                               BASIC EXTENSIONS
################################################################################
Base.length(s::CollectionSet)::Int = length(collection_sets(s))
Base.length(s::MultiDistributionSet)::Int = length(s.distribution)
Base.length(s::InfiniteScalarSet)::Int = 1

################################################################################
#                           SUPPORT VALIDITY METHODS
################################################################################
"""
    supports_in_set(supports::Union{Real, Vector{<:Real}, Array{<:Real, 2}},
                    set::AbstractInfiniteSet)::Bool

Used to check if `supports` are in the domain of `set`. Returns `true` if
`supports` are in domain of `set` and returns `false` otherwise.
This is primarily an internal method for performing checks but can be extended
for user-defined set types. Extending this is optional, but recommended where
possible. Note by fallback, this returns `true` for unrecognized set types such
that an error won't be thrown.
"""
function supports_in_set end

# IntervalSet
function supports_in_set(supports::Union{Real, Vector{<:Real}},
                         set::IntervalSet)::Bool
    min_support = minimum(supports)
    max_support = maximum(supports)
    if min_support < JuMP.lower_bound(set) || max_support > JuMP.upper_bound(set)
        return false
    end
    return true
end

# UnivariateDistribution set
function supports_in_set(supports::Union{Real, Vector{<:Real}},
                         set::UniDistributionSet)::Bool
    return all(Distributions.insupport(set.distribution, supports))
end

# MultivariateDistribution set
function supports_in_set(
    supports::Array{<:Real},
    set::MultiDistributionSet{<:Distributions.MultivariateDistribution{S}}
    )::Bool where {S <: Distributions.ValueSupport}
    if length(set.distribution) != size(supports, 1)
        error("Support dimensions does not match distribution dimensions.")
    end
    return all(Distributions.insupport(set.distribution, supports))
end

# CollectionSet
function supports_in_set(supports::Array{<:Real, 2}, set::CollectionSet)::Bool
    sets = collection_sets(set)
    if length(sets) != size(supports, 1)
        error("Support dimensions does not match CollectionSet dimensions.")
    end
    for i in eachindex(sets)
        if !supports_in_set(supports[i, :], sets[i])
            return false
        end
    end
    return true
end

# Fallback
function supports_in_set(supports, set::AbstractInfiniteSet)::Bool
    return true
end

################################################################################
#                              LOWER BOUND METHODS
################################################################################
"""
    JuMP.has_lower_bound(set::AbstractInfiniteSet)::Bool

Return `Bool` indicating if `set` has a lower bound that can be determined. This
should be extended for user-defined infinite sets. It defaults to `false`
for unrecognized set types.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> set = InfiniteSet(0, 1);

julia> has_lower_bound(set)
true
```
"""
function JuMP.has_lower_bound(set::AbstractInfiniteSet)::Bool # fallback
    return false
end

# IntervalSet
JuMP.has_lower_bound(set::IntervalSet)::Bool = true

# DistributionSet (Univariate)
JuMP.has_lower_bound(set::UniDistributionSet)::Bool = true

# CollectionSet
function JuMP.has_lower_bound(set::CollectionSet)::Bool
    for s in collection_sets(set)
        if !JuMP.has_lower_bound(s)
            return false
        end
    end
    return true
end

"""
    JuMP.lower_bound(set::AbstractInfiniteSet)::Union{Real, Vector{<:Real}}

Return the lower bound of `set` if one exists. This should be extended for
user-defined infinite sets if appropriate. Errors if `JuMP.has_lower_bound`
returns `false`. Extensions are enabled by `JuMP.has_lower_bound(set)` and
`JuMP.lower_bound(set)`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> set = InfiniteSet(0, 1);

julia> lower_bound(set)
0.0
```
"""
function JuMP.lower_bound(set::AbstractInfiniteSet) # fallback
    type = typeof(set)
    error("`JuMP.lower_bound` not defined for infinite set of type $type.")
end

# IntervalSet
JuMP.lower_bound(set::IntervalSet)::Float64 = set.lower_bound

# DistributionSet (Univariate)
function JuMP.lower_bound(set::UniDistributionSet)::Real
    return minimum(set.distribution)
end

# CollectionSet
function JuMP.lower_bound(set::CollectionSet)::Vector{<:Real}
    return [JuMP.lower_bound(i) for i in collection_sets(set)]
end

"""
    JuMP.set_lower_bound(set::AbstractInfiniteSet,
                         lower::Union{Real, Vector{<:Real}})::AbstractInfiniteSet

Set and return the lower bound of `set` if such an operation makes sense. Errors
if the type of `set` does not support this operation or has not been extended.
User-defined set types should extend this if appropriate.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> set = InfiniteSet(0, 1);

julia> set_lower_bound(set, 0.5)
[0.5, 1]
```
"""
function JuMP.set_lower_bound(set::AbstractInfiniteSet,
                              lower::Union{Real, Vector{<:Real}}) # fallback
    type = typeof(set)
    error("`JuMP.set_lower_bound` not defined for infinite set of type $type.")
end

# IntervalSet
function JuMP.set_lower_bound(set::IntervalSet, lower::Real)::IntervalSet
    return IntervalSet(lower, set.upper_bound)
end

# DistributionSet
function JuMP.set_lower_bound(set::Union{UniDistributionSet,
                                         MultiDistributionSet},
                              lower::Union{Real, Vector{<:Real}})
    error("Cannot set the lower bound of a distribution, try using " *
          "`Distributions.Truncated` instead.")
end

# CollectionSet
function JuMP.set_lower_bound(set::CollectionSet, lower::Vector{<:Real})::CollectionSet
    sets = collection_sets(set)
    new_sets = [JuMP.set_lower_bound(sets[i], lower[i]) for i in eachindex(sets)]
    return CollectionSet(new_sets)
end

################################################################################
#                              UPPER BOUND METHODS
################################################################################
"""
    JuMP.has_upper_bound(set::AbstractInfiniteSet)::Bool

Return `Bool` indicating if `set` has a upper bound that can be determined. This
should be extended for user-defined infinite sets. It defaults to `false`
for unrecognized set types.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> set = InfiniteSet(0, 1);

julia> has_upper_bound(set)
true
```
"""
function JuMP.has_upper_bound(set::AbstractInfiniteSet)::Bool # fallback
    return false
end

# IntervalSet
JuMP.has_upper_bound(set::IntervalSet)::Bool = true

# DistributionSet (Univariate)
JuMP.has_upper_bound(set::UniDistributionSet)::Bool = true

# CollectionSet
function JuMP.has_upper_bound(set::CollectionSet)::Bool
    for i in collection_sets(set)
        if !JuMP.has_upper_bound(i)
            return false
        end
    end
    return true
end

"""
    JuMP.upper_bound(set::AbstractInfiniteSet)::Union{Real, Vector{<:Real}}

Return the upper bound of `set` if one exists. This should be extended for
user-defined infinite sets if appropriate. Errors if `JuMP.has_upper_bound`
returns `false`. Extensions are enabled by `JuMP.has_upper_bound(set)` and
`JuMP.upper_bound(set)`.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> set = InfiniteSet(0, 1);

julia> upper_bound(set)
1.0
```
"""
function JuMP.upper_bound(set::AbstractInfiniteSet) # fallback
    type = typeof(set)
    error("`JuMP.upper_bound` not defined for infinite set of type $type.")
end

# IntervalSet
JuMP.upper_bound(set::IntervalSet)::Float64 = set.upper_bound

# DistributionSet (Univariate)
function JuMP.upper_bound(set::UniDistributionSet)::Real
    return maximum(set.distribution)
end

# CollectionSet
function JuMP.upper_bound(set::CollectionSet)::Vector{<:Real}
    return [JuMP.upper_bound(i) for i in collection_sets(set)]
end

"""
    JuMP.set_upper_bound(set::AbstractInfiniteSet,
                         upper::Real)::AbstractInfiniteSet

Set and return the upper bound of `set` if such an aoperation makes sense. Errors
if the type of `set` does not support this operation or has not been extended.
User-defined set types should extend this if appropriate.

**Example**
```jldoctest; setup = :(using InfiniteOpt, JuMP)
julia> set = InfiniteSet(0, 1);

julia> set_upper_bound(set, 0.5)
[0, 0.5]
```
"""
function JuMP.set_upper_bound(set::AbstractInfiniteSet, upper::Real) # fallback
    type = typeof(set)
    error("`JuMP.set_lower_bound` not defined for infinite set of type $type.")
end

# IntervalSet
function JuMP.set_upper_bound(set::IntervalSet, upper::Real)::IntervalSet
    return IntervalSet(set.lower_bound, upper)
end

# DistributionSet
function JuMP.set_upper_bound(set::Union{UniDistributionSet,
                                         MultiDistributionSet}, lower::Real)
    error("Cannot set the upper bound of a distribution, try using " *
          "`Distributions.Truncated` instead.")
end

# CollectionSet
function JuMP.set_upper_bound(set::CollectionSet, lower::Vector{<:Real})::CollectionSet
    sets = collection_sets(set)
    new_sets = [JuMP.set_upper_bound(sets[i], lower[i]) for i in eachindex(sets)]
    return CollectionSet(new_sets)
end

################################################################################
#                        SUPPORT AND LABEL GENERATION
################################################################################
"""
    AbstractSupportLabel

An abstract type for support label types. These are used to distinguish different
kinds of supports that are added to infinite parameters.
"""
abstract type AbstractSupportLabel end

"""
    All <: AbstractSupportLabel

This support label is unique in that it isn't associated with a particular set of
supports, but rather is used used to indicate that all supports should be used.
"""
struct All <: AbstractSupportLabel end

# Filler label for NoGenerativeSupports
struct _NoLabel <: AbstractSupportLabel end

"""
    PublicLabel <: AbstractSupportLabel

An abstract type used to denote that labels that should be given to the user by
default.
"""
abstract type PublicLabel <: AbstractSupportLabel end

"""
    UserDefined <: PublicLabel

A support label for supports that are supplied by the user directly to an infinite
parameter.
"""
struct UserDefined <: PublicLabel end

"""
    UniformGrid <: PublicLabel

A support label for supports that are generated uniformly accross a given interval.
"""
struct UniformGrid <: PublicLabel end

"""
    SampleLabel <: PublicLabel

An abstract type for labels of supports that are generated via some sampling technique.
"""
abstract type SampleLabel <: PublicLabel end

"""
    MCSample <: SampleLabel

A support label for supports that are generated via Monte Carlo Sampling.
"""
struct MCSample <: SampleLabel end

"""
    WeightedSample <: SampleLabel

A support label for supports that are generated by sampling from a statistical
distribution.
"""
struct WeightedSample <: SampleLabel end

"""
    Mixture <: PublicLabel

A support label for multi-dimensional supports that are generated from a variety
of methods.
"""
struct Mixture <: PublicLabel end

"""
    UniqueMeasure{S::Symbol} <: PublicLabel

A support label for supports that are provided from the `DiscreteMeasureData`
associated with a measure where a unique label is generated to distinguish those
supports. This is done by invoking [`generate_unique_label`](@ref).
"""
struct UniqueMeasure{S} <: PublicLabel end

"""
    InternalLabel <: AbstractSupportLabel

An abstract type for support labels that are associated with supports that should
not be reported to the user by default.
"""
abstract type InternalLabel <: AbstractSupportLabel end

"""
    generate_unique_label()::Type{UniqueMeasure}

Generate and return a unique support label for measures.
"""
function generate_unique_label()::DataType
    return UniqueMeasure{gensym()}
end

# Define default values of sig_digits and num_supports keywords
const DefaultSigDigits = 12
const DefaultNumSupports = 10

# a user interface of generate_support_values
"""
    generate_supports(set::AbstractInfiniteSet
                      [method::Type{<:AbstractSupportLabel}];
                      [num_supports::Int = DefaultNumSupports,
                      sig_digits::Int = DefaultSigDigits]
                      )::Tuple{Array{<:Real}, DataType}

Generate `num_supports` support values with `sig_digits` significant digits in
accordance with `set` and return them along with the correct generation label(s).
`IntervalSet`s generate supports uniformly with label `UniformGrid` and
distribution sets generate them randomly accordingly to the
underlying distribution. Moreover, `method` indicates the generation method that
should be used. These `methods` correspond to parameter support labels. Current
labels that can be used as generation methods include (but may not be defined
for certain set types):
- [`MCSample`](@ref): Uniformly distributed Monte Carlo samples.
- [`WeightedSample`](@ref): Monte Carlo samples that are weighted by an underlying PDF.
- [`UniformGrid`](@ref): Samples that are generated uniformly over the set domain.

Extensions that employ user-defined infinite set types and/or methods
should extend [`generate_support_values`](@ref) to enable this. Errors if the
`set` type and /or methods are unrecognized. This is intended as an internal
method to be used by methods such as [`generate_and_add_supports!`](@ref).
"""
function generate_supports(set::AbstractInfiniteSet;
                           num_supports::Int = DefaultNumSupports,
                           sig_digits::Int = DefaultSigDigits
                           )::Tuple
    return generate_support_values(set, num_supports = num_supports,
                                   sig_digits = sig_digits)
end

# 2 arguments
function generate_supports(set::AbstractInfiniteSet,
                           method::Type{<:AbstractSupportLabel};
                           num_supports::Int = DefaultNumSupports,
                           sig_digits::Int = DefaultSigDigits
                           )::Tuple
    return generate_support_values(set, method,
                                   num_supports = num_supports,
                                   sig_digits = sig_digits)
end

"""
    generate_support_values(set::AbstractInfiniteSet,
                            [method::Type{MyMethod} = MyMethod];
                            [num_supports::Int = DefaultNumSupports,
                            sig_digits::Int = DefaultSigDigits]
                            )::Tuple{Array{<:Real}, Symbol}

A multiple dispatch method for [`generate_supports`](@ref). This will return
a tuple where the first element are the supports and the second is their
label. This can be extended for user-defined infinite sets and/or generation
methods. When defining a new set type the default method dispatch should
make `method` an optional argument (making it the default). Otherwise, other
method dispatches for a given set must ensure that `method` is positional
argument without a default value (contrary to the definition above). Note that the
`method` must be a subtype of either [`PublicLabel`](@ref) or [`InternalLabel`](@ref).
"""
function generate_support_values(set::AbstractInfiniteSet,
                                 args...; kwargs...)
    if isempty(args)
        error("`generate_support_values` has not been extended for infinite sets " *
              "of type `$(typeof(set))`. This automatic support generation is not " *
              "implemented.")
    else
        error("`generate_support_values` has not been extended for infinite sets " *
              "of type `$(typeof(set))` with the generation method `$(args[1])`. " *
              "This automatic support generation is not implemented.")
    end
end

# IntervalSet and UniformGrid
function generate_support_values(set::IntervalSet,
                                 method::Type{UniformGrid} = UniformGrid;
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits,
                                 )::Tuple{Vector{<:Real}, DataType}
    lb = JuMP.lower_bound(set)
    ub = JuMP.upper_bound(set)
    new_supports = round.(range(lb, stop = ub, length = num_supports),
                          sigdigits = sig_digits)
    return new_supports, method
end

# IntervalSet and MCSample
function generate_support_values(set::IntervalSet,
                                 method::Type{MCSample};
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits,
                                 )::Tuple{Vector{<:Real}, DataType}
    lb = JuMP.lower_bound(set)
    ub = JuMP.upper_bound(set)
    dist = Distributions.Uniform(lb, ub)
    new_supports = round.(Distributions.rand(dist, num_supports),
                          sigdigits = sig_digits)
    return new_supports, method
end

# UniDistributionSet and MultiDistributionSet (with multivariate only)
function generate_support_values(
    set::Union{UniDistributionSet, MultiDistributionSet{<:Distributions.MultivariateDistribution}},
    method::Type{WeightedSample} = WeightedSample;
    num_supports::Int = DefaultNumSupports,
    sig_digits::Int = DefaultSigDigits
    )::Tuple{Array{<:Real}, DataType}
    dist = set.distribution
    new_supports = round.(Distributions.rand(dist, num_supports),
                          sigdigits = sig_digits)
    return new_supports, method
end

# UniDistributionSet and MCSample
function generate_support_values(
    set::UniDistributionSet,
    method::Type{MCSample};
    num_supports::Int = DefaultNumSupports,
    sig_digits::Int = DefaultSigDigits
    )::Tuple{Vector{Float64}, DataType}
    return generate_support_values(set, WeightedSample; num_supports = num_supports,
                                   sig_digits = sig_digits)[1], method # TODO use an unwieghted sample...
end

# MultiDistributionSet (matrix-variate distribution)
function generate_support_values(
    set::MultiDistributionSet{<:Distributions.MatrixDistribution},
    method::Type{WeightedSample} = WeightedSample;
    num_supports::Int = DefaultNumSupports,
    sig_digits::Int = DefaultSigDigits
    )::Tuple{Array{Float64, 2}, DataType}
    dist = set.distribution
    raw_supports = Distributions.rand(dist, num_supports)
    new_supports = Array{Float64}(undef, length(dist), num_supports)
    for i in 1:size(new_supports, 2)
        new_supports[:, i] = round.(reduce(vcat, raw_supports[i]),
                                    sigdigits = sig_digits)
    end
    return new_supports, method
end

# Generate the supports for a collection set
function _generate_collection_supports(set::CollectionSet, num_supports::Int,
                                       sig_digits::Int)::Array{Float64, 2}
    sets = collection_sets(set)
    # build the support array transpose to fill in column order (leverage locality)
    trans_supports = Array{Float64, 2}(undef, num_supports, length(sets))
    for i in eachindex(sets)
        @inbounds trans_supports[:, i] = generate_support_values(sets[i],
                                                   num_supports = num_supports,
                                                   sig_digits = sig_digits)[1]
    end
    return permutedims(trans_supports)
end

function _generate_collection_supports(set::CollectionSet,
                                       method::Type{<:AbstractSupportLabel},
                                       num_supports::Int,
                                       sig_digits::Int)::Array{Float64, 2}
    sets = collection_sets(set)
    # build the support array transpose to fill in column order (leverage locality)
    trans_supports = Array{Float64, 2}(undef, num_supports, length(sets))
    for i in eachindex(sets)
        @inbounds trans_supports[:, i] = generate_support_values(sets[i],
                                                   method,
                                                   num_supports = num_supports,
                                                   sig_digits = sig_digits)[1]
    end
    return permutedims(trans_supports)
end

# CollectionSet (IntervalSets)
function generate_support_values(set::CollectionSet{IntervalSet},
                                 method::Type{UniformGrid} = UniformGrid;
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(set, num_supports, sig_digits)
    return new_supports, method
end

function generate_support_values(set::CollectionSet{IntervalSet},
                                 method::Type{MCSample};
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(set, method, num_supports, sig_digits)
    return new_supports, method
end

# CollectionSet (UniDistributionSets)
function generate_support_values(set::CollectionSet{<:UniDistributionSet},
                                 method::Type{WeightedSample} = WeightedSample;
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(set, num_supports, sig_digits)
    return new_supports, method
end

# CollectionSet (InfiniteScalarSets)
function generate_support_values(set::CollectionSet,
                                 method::Type{Mixture} = Mixture;
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(set, num_supports, sig_digits)
    return new_supports, method
end

# CollectionSet (InfiniteScalarSets) using purely MC sampling
# this is useful for measure support generation
function generate_support_values(set::CollectionSet,
                                 method::Type{MCSample};
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits
                                 )::Tuple{Array{<:Real}, DataType}
    new_supports = _generate_collection_supports(set, method, num_supports, sig_digits)
    return new_supports, method
end

# For label All: dispatch to default methods
function generate_support_values(set::AbstractInfiniteSet, ::Type{All};
                                 num_supports::Int = DefaultNumSupports,
                                 sig_digits::Int = DefaultSigDigits)
    return generate_support_values(set, num_supports = num_supports,
                                   sig_digits = sig_digits)
end


text:
using Gaston
g = Gaston;

set_terminal("x11");

c = g.CurveConf();
c.plotstyle = "pm3d";

x = -15:0.33:15;
y = -15:0.33:15;
Z = g.meshgrid(x, y, (x,y) -> sin(sqrt(x .* x+y .* y))
					 / sqrt(x .* x+y .* y));
g.addcoords(x,y,Z,c);

a = g.AxesConf();
a.title = "3D: Sombrero";
g.addconf(a);
g.llplot()


text:
using Crayons.Box

println(WHITE_FG, BLACK_BG, "Normal")
println(WHITE_BG, BLACK_FG, "Reversed")
println(WHITE_FG, BLACK_BG, "Normal")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#Broadcasting machinery
#######################

# YOU SHALL NOT PASS:
# For real though, some of this is far from pretty.

valunwrap(::Val{x}) where x = x
#This is directly from the Broadcasting interface docs
struct GBVectorStyle <: Broadcast.AbstractArrayStyle{1} end
struct GBMatrixStyle <: Broadcast.AbstractArrayStyle{2} end
Base.BroadcastStyle(::Type{<:GBVector}) = GBVectorStyle()
Base.BroadcastStyle(::Type{<:GBMatrix}) = GBMatrixStyle()
Base.BroadcastStyle(::Type{<:Transpose{T, <:GBMatrix} where T}) = GBMatrixStyle()
#
GBVectorStyle(::Val{0}) = GBVectorStyle()
GBVectorStyle(::Val{1}) = GBVectorStyle()
GBVectorStyle(::Val{2}) = GBMatrixStyle()
GBVectorStyle(::Val{N}) where N = Broadcast.DefaultArrayStyle{N}()
GBMatrixStyle(::Val{0}) = GBMatrixStyle()
GBMatrixStyle(::Val{1}) = GBMatrixStyle()
GBMatrixStyle(::Val{2}) = GBMatrixStyle()
GBMatrixStyle(::Val{N}) where N = Broadcast.DefaultArrayStyle{N}()
Broadcast.BroadcastStyle(::GBMatrixStyle, ::GBVectorStyle) = GBMatrixStyle()

function Base.similar(
    bc::Broadcast.Broadcasted{GBMatrixStyle},
    ::Type{ElType}
) where {ElType}
    return GBMatrix{ElType}(axes(bc))
end

function Base.similar(
    bc::Broadcast.Broadcasted{GBVectorStyle},
    ::Type{ElType}
) where {ElType}
    return GBVector{ElType}(axes(bc))
end

#Find the modifying version of a function.
modifying(::typeof(mul)) = mul!
modifying(::typeof(eadd)) = eadd!
modifying(::typeof(emul)) = emul!

@inline function Base.copy(bc::Broadcast.Broadcasted{GBMatrixStyle})
    f = bc.f
    l = length(bc.args)
    if l == 1
        x = first(bc.args)
        if x isa Broadcast.Broadcasted
            x = copy(x)
        end
        return map(f, x)
    else
        left = first(bc.args)
        right = last(bc.args)
        if left isa Base.RefValue{typeof(^)}
            f = ^
            left = bc.args[2]
            right = valunwrap(right[])
        end
        if left isa Broadcast.Broadcasted
            left = copy(left)
        end
        if right isa Broadcast.Broadcasted
            right = copy(right)
        end
        if left isa GBArray && right isa GBArray
            add = defaultadd(f)
            return add(left, right, f)
        else
            return map(f, left, right)
        end
    end
end
mutatingop(::typeof(emul)) = emul!
mutatingop(::typeof(eadd)) = eadd!
mutatingop(::typeof(map)) = map!
@inline function Base.copyto!(C::GBArray, bc::Broadcast.Broadcasted{GBMatrixStyle})
    l = length(bc.args)
    if l == 1
        x = first(bc.args)
        if bc.f === Base.identity
            C[:,:, accum=second] = x
            return C
        end
        return map!(bc.f, C, x; accum=second)
    else

        left = first(bc.args)
        right = last(bc.args)
        # handle annoyances with the pow operator
        if left isa Base.RefValue{typeof(^)}
            f = ^
            left = bc.args[2]
            right = valunwrap(right[])
        end
        # TODO: This if statement should probably be *inside* one of the inner ones to avoid duplication.
        if left === C
            if !(right isa Broadcast.Broadcasted)
                # This should be something of the form A .<op>= <expr> or A .= A .<op> <expr> which are equivalent.
                # this will be done by a subassign
                C[:,:, accum=bc.f] = right
                return C
            else
                # The form A .<op>= expr
                # but not of the form A .= C ... B.
                accum = bc.f
                f = right.f
                if length(right.args) == 1
                    # Should be catching expressions of the form A .<op>= <op>.(B)
                    subarg = first(right.args)
                    if subarg isa Broadcast.Broadcasted
                        subarg = copy(subarg)
                    end
                    return map!(f, C, subarg; accum)
                else
                    # Otherwise we know there's two operands on the LHS so we have A .<op>= C .<op> B
                    # Or a generalization with any compound *lazy* RHS.
                    (subargleft, subargright) = right.args
                    # subargleft and subargright are C and B respectively.
                    # If they're further nested broadcasts we can't fuse them, so just copy.
                    subargleft isa Broadcast.Broadcasted && (subargleft = copy(subargleft))
                    subargright isa Broadcast.Broadcasted && (subargright = copy(subargright))
                    if subargleft isa GBArray && subargright isa GBArray
                        add = mutatingop(defaultadd(f))
                        return add(C, subargleft, subargright, f; accum)
                    else
                        return map!(f, C, subargleft, subargright; accum)
                    end
                end
            end
        else
            # Some expression of the form A .= C .<op> B or a generalization
            # excluding A .= A .<op> <expr>, since that is captured above.
            if left isa Broadcast.Broadcasted
                left = copy(left)
            end
            if right isa Broadcast.Broadcasted
                right = copy(right)
            end
            if left isa GBArray && right isa GBArray
                add = mutatingop(defaultadd(f))
                return add(C, left, right, f)
            else
                return map!(C, f, left, right; accum=second)
            end
        end
    end
end

@inline function Base.copy(bc::Broadcast.Broadcasted{GBVectorStyle})
    f = bc.f
    l = length(bc.args)
    if l == 1
        x = first(bc.args)
        if x isa Broadcast.Broadcasted
            x = copy(x)
        end
        return map(f, x)
    else
        left = first(bc.args)
        right = last(bc.args)
        if left isa Base.RefValue{typeof(^)}
            f = ^
            left = bc.args[2]
            right = valunwrap(right[])
        end
        if left isa Broadcast.Broadcasted
            left = copy(left)
        end
        if right isa Broadcast.Broadcasted
            right = copy(right)
        end
        if left isa GBArray && right isa GBArray
            add = defaultadd(f)
            return add(left, right, f)
        else
            return map(f, left, right)
        end
    end
end

#Elementwise Broadcasts
#####################
#function Base.broadcasted(::typeof(*), A::GBArray, B::GBArray)
#    emul(A, B, BinaryOps.TIMES)
#end
#
#function Base.broadcasted(::typeof(+), A::GBArray, B::GBArray)
#    eadd(A, B, BinaryOps.PLUS)
#end
#
#function Base.broadcasted(::typeof(-), A::GBArray, B::GBArray)
#    eadd(A, B, BinaryOps.MINUS)
#end
#
##TODO: fix tricky gotchas, this will do type-specific (ie sometimes integer) division.
#function Base.broadcasted(::typeof(/), A::GBArray, B::GBArray)
#    emul(A, B, BinaryOps.DIV)
#end
#
##TODO: fix tricky gotchas, this will do type-specific (ie sometimes integer) division.
#function Base.broadcasted(::typeof(\), A::GBArray, B::GBArray)
#    emul(A, B, BinaryOps.RDIV)
#end
#
#function Base.broadcasted(::typeof(==), A::GBArray, B::GBArray)
#    eadd(A, B, BinaryOps.EQ)
#end
#
#function Base.broadcasted(::typeof(!=), A::GBArray, B::GBArray)
#    eadd(A, B, BinaryOps.NE)
#end
#
#function Base.broadcasted(::typeof(<), A::GBArray, B::GBArray)
#    eadd(A, B, BinaryOps.LT)
#end
#
#function Base.broadcasted(::typeof(>), A::GBArray, B::GBArray)
#    eadd(A, B, BinaryOps.GT)
#end
#
#function Base.broadcasted(::typeof(<=), A::GBArray, B::GBArray)
#    eadd(A, B, BinaryOps.LE)
#end
#
#function Base.broadcasted(::typeof(>=), A::GBArray, B::GBArray)
#    eadd(A, B, BinaryOps.GE)
#end
#
#function Base.broadcasted(::typeof(^), A::GBArray, B::GBArray)
#    emul(A, B, BinaryOps.POW)
#end
#
##Map Broadcasts
########################
#function Base.broadcasted(::typeof(*), u::GBArray, x::valid_union)
#    map(BinaryOps.TIMES, u, x)
#end
#function Base.broadcasted(::typeof(*), x::valid_union, u::GBArray)
#    map(BinaryOps.TIMES, x, u)
#end
#
#function Base.broadcasted(::typeof(/), u::GBArray, x::valid_union)
#    map(BinaryOps.DIV, u, x)
#end
#function Base.broadcasted(::typeof(/), x::valid_union, u::GBArray)
#    map(BinaryOps.DIV, x, u;)
#end
#
#function Base.broadcasted(::typeof(^), u::GBArray, x::valid_union)
#    map(BinaryOps.POW, u, x)
#end
#function Base.broadcasted(::typeof(^), x::valid_union, u::GBArray)
#    map(BinaryOps.POW, x, u)
#end
#
#
#function Base.broadcasted(::typeof(-), u::GBArray, x::valid_union)
#    map(BinaryOps.MINUS, u, x)
#end
#function Base.broadcasted(::typeof(-), x::valid_union, u::GBArray)
#    map(BinaryOps.MINUS, x, u)
#end
#function Base.broadcasted(::typeof(+), u::GBArray, x::valid_union)
#    map(BinaryOps.PLUS, u, x)
#end
#function Base.broadcasted(::typeof(+), x::valid_union, u::GBArray)
#    map(BinaryOps.PLUS, x, u)
#end


text:
"""
    get_bea_parameterlist(b::BeaData.Bea, dataset::String)

Return a `Dict` of parameters for `dataset`.

Arguments
---------
* `b` -- a [`Bea`](@ref) connection
* `dataset` -- String indicating the dataset ID.

"""
function get_bea_parameterlist(b::Bea, dataset::String)
    url = b.url
    key = b.key
    bea_method = "GetParameterList"

    querydict = Dict("UserID" => key,
                     "Method" => bea_method,
                     "DatasetName" => dataset,
                     "ResultFormat" => "JSON")

    response_json = bea_query(url, querydict)
    response_dict = response_json["BEAAPI"]["Results"]["Parameter"]

    return response_dict
end