text:
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, HydroDispatchRunOfRiverOutage},
    ::Type{S},
) where {H <: PSY.HydroGen, S <: PM.AbstractPowerModel}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerVariable,
        devices,
        HydroDispatchRunOfRiverOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ReactivePowerVariable,
        devices,
        HydroDispatchRunOfRiverOutage(),
    )

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ReactivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end

"""
Construct model for HydroGen with RunOfRiver Dispatch Formulation
with only Active Power.
"""
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, HydroDispatchRunOfRiverOutage},
    ::Type{S},
) where {H <: PSY.HydroGen, S <: PM.AbstractActivePowerModel}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerVariable,
        devices,
        HydroDispatchRunOfRiverOutage(),
    )

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end

function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, HydroDispatchReservoirStorageOutage},
    ::Type{S},
) where {H <: PSY.HydroEnergyReservoir, S <: PM.AbstractPowerModel}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ReactivePowerVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.SpillageVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyShortageVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergySurplusVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )

    # Initial Conditions
    PSI.storage_energy_initial_condition!(
        optimization_container,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    # Energy Balance Constraint
    PSI.add_constraints!(
        optimization_container,
        PSI.EnergyBalanceConstraint,
        PSI.EnergyVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.energy_target_constraint!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end

"""
Construct model for HydroGen with ReservoirStorage Dispatch Formulation
with only Active Power
"""
function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, HydroDispatchReservoirStorageOutage},
    ::Type{S},
) where {H <: PSY.HydroEnergyReservoir, S <: PM.AbstractActivePowerModel}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    # Variables

    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.SpillageVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyShortageVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergySurplusVariable,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )

    # Initial Conditions
    PSI.storage_energy_initial_condition!(
        optimization_container,
        devices,
        HydroDispatchReservoirStorageOutage(),
    )
    # Energy Balance Constraint
    PSI.add_constraints!(
        optimization_container,
        PSI.EnergyBalanceConstraint,
        PSI.EnergyVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.energy_target_constraint!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end

function PSI.construct_device!(
    optimization_container::PSI.OptimizationContainer,
    sys::PSY.System,
    model::PSI.DeviceModel{H, HydroDispatchPumpedStoragewReservationOutage},
    ::Type{S},
) where {H <: PSY.HydroPumpedStorage, S <: PM.AbstractActivePowerModel}
    devices = PSI.get_available_components(H, sys)

    if !PSI.validate_available_devices(H, devices)
        return
    end

    # Variables
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerInVariable,
        devices,
        HydroDispatchPumpedStoragewReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ActivePowerOutVariable,
        devices,
        HydroDispatchPumpedStoragewReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyVariableUp,
        devices,
        HydroDispatchPumpedStoragewReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.EnergyVariableDown,
        devices,
        HydroDispatchPumpedStoragewReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.SpillageVariable,
        devices,
        HydroDispatchPumpedStoragewReservationOutage(),
    )
    PSI.add_variables!(
        optimization_container,
        PSI.ReserveVariable,
        devices,
        HydroDispatchPumpedStoragewReservationOutage(),
    )

    # Constraints
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerOutVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.RangeConstraint,
        PSI.ActivePowerInVariable,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )

    # Initial Conditions
    PSI.storage_energy_initial_condition!(
        optimization_container,
        devices,
        HydroDispatchPumpedStoragewReservationOutage(),
    )

    # Energy Balanace limits
    PSI.add_constraints!(
        optimization_container,
        PSI.EnergyBalanceConstraint,
        PSI.EnergyVariableUp,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.add_constraints!(
        optimization_container,
        PSI.EnergyBalanceConstraint,
        PSI.EnergyVariableDown,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    outage_constraints!(
        optimization_container,
        devices,
        model,
        S,
        PSI.get_feedforward(model),
    )
    PSI.feedforward!(optimization_container, devices, model, PSI.get_feedforward(model))

    # Cost Function
    PSI.cost_function!(optimization_container, devices, model, S, nothing)

    return
end


text:
"""
    biconnectedComponents(network, ignoreTrivial=false)

Calculate biconnected components (aka "blobs") using Tarjan's algorithm:
the output is an array of arrays of edges.
These blobs are returned in post-order, but within a blob,
edges are *not* necessarily sorted in topological order.
If `ignoreTrivial` is true, trivial components (of a single edge)
are not returned.
The network is assumed to be connected.

**Warnings**: for nodes, fields `k`, `inCycle`, and `prev`
are modified during the algorithm. They are used to store the
node's "index" (time of visitation), "lowpoint", and the node's
"parent", as defined by the order in which nodes are visited.

References:
- p. 153 of [Tarjan (1972)](http://epubs.siam.org/doi/pdf/10.1137/0201010).
  Depth first search and linear graph algorithms,
  SIAM Journal on Computing, 1(2):146-160
- on [geeksforgeeks](https://www.geeksforgeeks.org/biconnected-components/),
  there is an error (as of 2018-01-30):
  `elif v != parent[u] and low[u] > disc[v]:` (python version)
  should be replaced by
  `elif v != parent[u] and disc[u] > disc[v]:`
- nice explanation at this
  [url](https://www.cs.cmu.edu/~avrim/451f12/lectures/biconnected.pdf)
"""
function biconnectedComponents(net, ignoreTrivial=false::Bool)
    for n in net.node
        n.inCycle = -1 # -1 for missing, node not visited yet
        n.k = -1       # inCycle = lowpoint
        n.prev = nothing # parent: will be a node later
    end
    S = Edge[] # temporary stack
    blobs = Vector{Edge}[] # will contain the blobs
    biconnectedComponents(net.node[net.root], [0], S, blobs, ignoreTrivial)
    # if stack not empty: create last connected component
    if length(S)>0
        #println("stack not empty at the end.")
        bb = Edge[]
        while length(S)>0
            e2 = pop!(S)
            push!(bb, e2)
        end
        if !ignoreTrivial || length(bb)>1
            push!(blobs, bb)
        end
    end
    return blobs
end

"""
    biconnectedComponents(node, index, S, blobs, ignoreTrivial)

Helper recursive function starting at a node (not a network).
`index` is an array containing a single integer, thus mutable:
order in which nodes are visited.
"""
function biconnectedComponents(node, index, S, blobs, ignoreTrivial)
    #println("\nentering biconnect, index=$(index[1])")
    children = 0
    # set depth index for v to the smallest unused index
    node.k = index[1]       # index / disc
    node.inCycle = index[1] # lowpoint / lowlink
    #println(" for node $(node.number): k=$(node.k), low=$(node.inCycle)")
    index[1] += 1
    #print(" stack: "); @show [e.number for e in S]

    for e in node.edge # each (v, w) in E do
        w = (e.node[1] == node) ? e.node[2] : e.node[1]
        #println(" w: $(w.number) along edge $(e.number)")
        if w.k == -1 # w not yet visited
            w.prev = node # set parent of w to current node
            #println(" w's parent = $(w.prev.number)")
            children += 1
            push!(S, e)
            #println(" edge $(e.number) on stack, intermediate step for node=$(node.number)")
            biconnectedComponents(w, index, S, blobs, ignoreTrivial)
            # check if subtree rooted at w has connection to
            # one of the ancestors of node
            # Case 1 -- per Strongly Connected Components Article
            node.inCycle = min(node.inCycle, w.inCycle) # lowpoint
            #print(" case 1, low=$(node.inCycle) for node $(node.number); ")
            #println("low=$(w.inCycle) for w $(w.number)")
            # if node is an articulation point: pop until e
            if (node.prev === nothing && children > 1) ||
               (node.prev !== nothing && w.inCycle >= node.k)
                # node is either root or an articulation point
                # start a new strongly connected component
                bb = Edge[]
                while length(S)>0
                    e2 = pop!(S)
                    #println(" popping: edge $(e2.number)")
                    push!(bb, e2)
                    e2 != e || break
                end
                if !ignoreTrivial || length(bb)>1
                    push!(blobs, bb)
                end
            end
        elseif w != node.prev && node.k > w.k
            # e is a back edge, not cross edge. Case 2 in article.
            node.inCycle = min(node.inCycle, w.k)
            #println(" case 2, node $(node.number): low=$(node.inCycle)")
            push!(S, e)
            #println(" edge $(e.number) now on stack, final step for node $(node.number)")
        else
            #println(" nothing to do from node $(node.number) to w $(w.number) along edge $(e.number)")
        end
    end
    return nothing
end

"""
    blobInfo(network, ignoreTrivial=true)

Calculate the biconnected components (blobs) using function
[`biconnectedComponents`](@ref) then:
- set node field `isExtBadTriangle` to true at the root of each
  non-trivial blob (and at the network root), false otherwise.
  (a better name for the field would be something like "isBlobRoot".)
- output:
  1. array of nodes that are the roots of each non-trivial blob,
     and the network root. If the root of the full network is
     not part of a non-trivial blob, a corresponding blob is
     added to the list.
  2. array of arrays: for each non-trivial blob,
     array of major hybrid edges in that blob.
  3. array of arrays: same as #2 but for minor hybrid edges,
     with hybrids listed in the same order, for each blob.

Blobs are ordered in reverse topological ordering
(aka post order).
If `ignoreTrivial` is true, trivial components are ignored.

keyword argument: `checkPreorder`, true by default. If false,
the `isChild1` edge field and the `net.nodes_changed` network field
are supposed to be correct.
"""
function blobInfo(net, ignoreTrivial=true::Bool;
    checkPreorder=true::Bool)
    if checkPreorder
      directEdges!(net) # update isChild1, needed for preorder
      preorder!(net) # creates / updates net.nodes_changed
    end
    bcc = biconnectedComponents(net, ignoreTrivial)
    bccRoots = Node[]         # one root node for each blob
    bccMajor = Vector{Edge}[] # one array for each blob
    bccMinor = Vector{Edge}[]
    for bicomp in bcc
        bccMa = Edge[]
        jmin = length(net.node)
        for edge in bicomp
            n = getParent(edge)
            j = something(findfirst(x -> x===n, net.nodes_changed), 0)
            jmin = min(j, jmin)
            if edge.hybrid && edge.isMajor
                push!(bccMa, edge)
            end
        end
        push!(bccRoots, net.nodes_changed[jmin])
        push!(bccMajor, bccMa)
        bccmi = Edge[] # find minor hybrid edges, in same order
        for edge in bccMa
            e = getPartner(edge)
            !e.isMajor || @warn "major edge $(edge.number) has a major partner: edge $(e.number)"
            push!(bccmi, e)
        end
        push!(bccMinor, bccmi)
    end
    # add the network root, if it was in a trivial bi-component (no hybrids)
    # blobs in post-order, so if there the root is there, it's the last blob
    if length(bcc)==0 || bccRoots[end]!=net.node[net.root]
        push!(bccRoots,net.node[net.root])
        push!(bccMajor, Edge[])
        push!(bccMinor, Edge[])
    end
    # update `isExtBadTriangle` to mark nodes that are blob roots:
    # these nodes will serve as dummy leaves when reached from other blobs
    for n in net.node n.isExtBadTriangle = false; end
    for r in bccRoots r.isExtBadTriangle = true;  end
    return(bccRoots, bccMajor, bccMinor)
end

"""
    blobDecomposition!(network)
    blobDecomposition(network)

Find blobs using [`biconnectedComponents`](@ref); find their roots
using [`blobInfo`](@ref); create a forest in the form of a
disconnected network (for efficiency), by deconnecting the
root of each non-trivial blob from its parent.
The root of each blob corresponds to a new leaf
(in another tree of the forest):
the number of the blob's root is given to the newly created leaf.

The first (bang) version modifies the network and returns
the array of blob roots. The second version copies the network
then returns a tuple: the forest and the blob roots.
"""
function blobDecomposition(net)
    net2 = deepcopy(net)
    blobR = blobDecomposition!(net2)
    return net2, blobR
end
function blobDecomposition!(net)
    nextnumber = maximum([n.number for n in net.node])+1
    blobR, tmp, tmp = blobInfo(net, true) # true: ignore trivial single-edge blobs
    for r in blobR
        for e in r.edge
            r == e.node[e.isChild1 ? 1 : 2] || continue
            removeEdge!(r,e) # detach edge e from root r
            removeNode!(r,e) # detach root r from edge e
            dummyleaf = Node(nextnumber, true) # true: leaf
            nextnumber += 1
            dummyleaf.name = string("dummy ", r.number)
            setEdge!(dummyleaf, e) # attach e to new dummy leaf
            setNode!(e,dummyleaf)  # attach new leaf to edge e
            e.isChild1 = false
            pushNode!(net, dummyleaf)
            break
        end
    end
    return blobR
end

"""
    treeedgecomponents(net::HybridNetwork)

Return the tree-edge components of the semidirected network as a `membership`
dictionary `Node => Int`. Nodes with the same membership integer
value are in the same tree-edge component.
The tree-edge components of a network are the connected components of
the network when all hybrid edges are removed.

A `RootMismatch` error is thrown if there exists a cycle in any of the tree-edge
components, or if a tree-edge component has more than one "entry" hybrid node.

Warnings:
- since `Node`s are mutable, the network should not be modified
  until usage of the output `membership` dictionary is over.
- the component IDs are not predicable, but will be consecutive integers
  from 1 to the number of components.
"""
function treeedgecomponents(net::HybridNetwork)
    # partition nodes into tree-edge components (TECs)
    nodes = net.node
    n = length(nodes)
    unvisited = Set(nodes)
    dfs_stack = Vector{Node}()  # stack for iterative depth-first search over tree edges
    dfs_parent = Dict{Node, Node}() # dfs_parent[node] = node's parent in DFS tree
    # membership[node] = id of undirected component that the node belongs to
    membership = Dict{Node, Int}()
    cur_id = 0                  # undirected component id

    while !isempty(unvisited)
        # loop over undirected components
        # start with an unvisited node, DFS with undirected edges
        cur_id += 1
        node = pop!(unvisited) # unpredictable ordering: because unvisited is a set
        push!(dfs_stack, node)
        curnode = node
        dfs_parent[curnode] = curnode
        entrynode = nothing
        while !isempty(dfs_stack)
            # DFS loop over one tree-edge component
            curnode = pop!(dfs_stack)
            delete!(unvisited, curnode)
            membership[curnode] = cur_id
            for e in curnode.edge
                if !e.hybrid
                    # for tree edge, do DFS, check for undirected cycles
                    nextnode = getOtherNode(e, curnode)
                    # if run into visited node (other than parent), then component has cycle
                    if nextnode !== dfs_parent[curnode]
                        if !in(nextnode, unvisited)
                            throw(RootMismatch(
                                "Undirected cycle exists, starting at node number $(nextnode.number)"))
                        else
                            dfs_parent[nextnode] = curnode
                            push!(dfs_stack, nextnode)
                        end
                    end
                else # for hybrid edge, check there is at most one entry node into the TEC
                    if curnode === getChild(e)
                        if isnothing(entrynode)
                            entrynode = curnode
                        elseif entrynode !== curnode
                            throw(RootMismatch(
                                """Multiple entry nodes for tree-edge component, numbered:
                                $(entrynode.number) and $(curnode.number)"""))
                        end
                    end
                end
            end
        end
    end

    return membership
end

"""
    checkroot!(net)
    checkroot!(net::HybridNetwork, membership::Dict{Node, Int})

Set the root of `net` to an appropriate node and update the edges `containRoot`
field appropriately, using the `membership` output by [`treeedgecomponents`](@ref).
A node is appropriate to serve as root if it belongs in the
root tree-edge component, that is, the root of the tree-edge component graph.

- If the current root is appropriate, it is left as is. The direction of
  edges (via `isChild1`) is also left as is, assuming it was in synch with
  the existing root.
- Otherwise, the root is set to the first appropriate node in `net.node`,
  that is not a leaf. Then edges are directed away from this root.

A `RootMismatch` error is thrown if `net` is not a valid semidirected
phylogenetic network (i.e. it is not possible to root the network in a way
compatible with the given hybrid edges).

Output: the `membership` ID of the root component.
The full set of nodes in the root component can be obtained as shown below.
Warning: only use the output component ID after calling the second version
`checkroot!(net, membership)`.

```jldoctest
julia> net = readTopology("(#H1:::0.1,#H2:::0.2,(((b)#H1)#H2,a));");

julia> membership = treeedgecomponents(net);

julia> rootcompID = checkroot!(net, membership);

julia> rootcomp = keys(filter(p -> p.second == rootcompID, membership));

julia> sort([n.number for n in rootcomp]) # number of nodes in the root component
3-element Vector{Int64}:
 -3
 -2
  4
```
"""
function checkroot!(net::HybridNetwork)
    membership = treeedgecomponents(net) # dict Node => Int
    return checkroot!(net, membership)
end
function checkroot!(net::HybridNetwork, membership::Dict{Node,Int})
    # do not modify nodes or edges until we are done with membership
    nodes = net.node
    ncomp = maximum(values(membership)) # TECs are numbered 1,2,...,ncomp
    #= 1. construct TEC graph: directed graph in which
    vertices = TECs of original network, and
    edge Ci --> Cj for each original hybrid edge: node in Ci --> node in Cj.
    =#
    tecG = [Set{Int}() for comp in 1:ncomp] # tecG[i] = set of children of ith TEC
    noparent = trues(ncomp) # will stay true for TECs with no parent
    for e in net.edge
        e.hybrid || continue # skip tree edges
        up, down = e.isChild1 ? (e.node[2], e.node[1]) : (e.node[1], e.node[2])
        uc_up, uc_down = membership[up], membership[down]
        noparent[uc_down] = false
        push!(tecG[uc_up], uc_down)
    end
    # 2. check that the TEC graph has a single root
    tec_root = findfirst(noparent)
    if isnothing(tec_root) # all TECs have a parent: noparent are all false
        throw(RootMismatch(
            "Semidirected cycle exists: no component has in-degree 0"))
    elseif sum(noparent) > 1
        r1 = findfirst(noparent)
        r2 = findlast(noparent)
        nr1 = findfirst(n -> membership[n] == r1, nodes)
        nr2 = findfirst(n -> membership[n] == r2, nodes)
        throw(RootMismatch("nodes number $(nodes[nr1].number) and $(nodes[nr2].number) have no common ancestor"))
    end

    # 3. topological sort: check the TEC graph has no cycles, that is, is a DAG
    indeg = zeros(Int, ncomp)  # in-degree of vertex in TEC graph
    for uc in tecG
        for i in uc
            indeg[i] += 1
        end
    end
    headstack = [tec_root] # stack of TEC vertices of degree 0 in topological sort
    while !isempty(headstack)
        uc = pop!(headstack)
        indeg[uc] -= 1
        for i in tecG[uc]
            indeg[i] -= 1
            if indeg[i] == 0
                push!(headstack, i)
            end
        end
    end
    cyclehead = findfirst(indeg .!= -1)
    isnothing(cyclehead) || throw(RootMismatch(
        """Semidirected cycle exists, starting at TEC containing node number $(nodes[findfirst(n -> membership[n] == cyclehead, nodes)].number)"""))

    # 4. original network: reset the network root if needed,
    #    and update the edges' containRoot accordingly.
    # NOT done: build the root component to return it. Instead: return tec_root
    #   Set(node for node = nodes if membership[node] == tec_root)
    #   Set(node for (node,comp) in membership if comp == tec_root)
    #rootcomp = keys(filter(p -> p.second == tec_root, membership))
    curroot = nodes[net.root]
    if membership[curroot] == tec_root
        # update containRoot only: true for edges in or out of the root TEC
        for e in net.edge
            e.containRoot = (membership[getParent(e)] == tec_root)
        end
    else
        net.root = findfirst(n -> (!n.leaf && membership[n] == tec_root), nodes)
        directEdges!(net) # also updates containRoot of all edges
    end
    return tec_root # return rootcomp
end

ext: jl


text:
using BinaryBuilder

name = "Shoco"
version = v"2015.10.8"  # No tagged versions, this is the date of the commit used

sources = [
    GitSource("https://github.com/Ed-von-Schleck/shoco.git",
              "4dee0fc850cdec2bdb911093fe0a6a56e3623b71"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/shoco/
mkdir -p ${libdir}
${CC} shoco.c -o "${libdir}/libshoco.${dlext}" -fPIC -std=c99 -shared
"""

platforms = supported_platforms(; experimental=true)

products = [
    LibraryProduct("libshoco", :libshoco),
]

dependencies = []

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
               julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function get_snake(;
    timestep=0.01,
    gravity=[0.0; 0.0; -9.81],
    friction_coefficient=0.8,
    contact=true,
    contact_type=:nonlinear,
    spring=0.0,
    damper=0.0,
    num_bodies=2,
    joint_type=:Spherical,
    length=1.0,
    radius=0.05,
    T=Float64)

    # Parameters
    ex = [1.0; 0.0; 0.0]
    ey = [0.0; 1.0; 0.0]
    ez = [0.0; 0.0; 1.0]

    vert11 = [0.0; 0.0; length / 2.0]
    vert12 = -vert11

    # Links
    origin = Origin{T}()
    bodies = [Box(3 * radius, 2 * radius, length, length, color=RGBA(0.9, 0.9, 0.9)) for i = 1:num_bodies]

    # Constraints
    jointb1 = JointConstraint(Floating(origin, bodies[1],
        spring=0.0,
        damper=0.0))

    if num_bodies > 1
        joints = [JointConstraint(Prototype(joint_type, bodies[i - 1], bodies[i], ex;
            parent_vertex=vert12,
            child_vertex=vert11,
            spring=spring,
            damper=damper)) for i = 2:num_bodies]
        joints = [jointb1; joints]
    else
        joints = [jointb1]
    end

    if contact
        n = num_bodies
        normal = [[0.0; 0.0; 1.0] for i = 1:n]
        friction_coefficient = friction_coefficient * ones(n)

        contacts1 = contact_constraint(bodies, normal,
            friction_coefficient=friction_coefficient,
            contact_origins=fill(vert11, n),
            contact_type=contact_type) # we need to duplicate point for prismatic joint for instance
        contacts2 = contact_constraint(bodies, normal,
            friction_coefficient=friction_coefficient,
            contact_origins=fill(vert12, n),
            contact_type=contact_type)
        mech = Mechanism(origin, bodies, joints, [contacts1; contacts2],
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    else
        mech = Mechanism(origin, bodies, joints,
            gravity=gravity,
            timestep=timestep,
            spring=spring,
            damper=damper)
    end
    return mech
end

function initialize_snake!(mechanism::Mechanism{T};
    base_position=[0.0, -0.5, 0.0],
    base_orientation=RotX(0.6 * π),
    base_linear_velocity=zeros(3),
    base_angular_velocity=zeros(3),
    relative_linear_velocity=zeros(3),
    relative_angular_velocity=zeros(3)) where T

    pbody = mechanism.bodies[1]
    h = pbody.shape.xyz[3]
    vert11 = [0.0; 0.0; h / 2.0]
    vert12 = -vert11

    # set position and velocities
    set_maximal_configurations!(mechanism.origin, pbody,
        child_vertex=base_position,
        Δq=base_orientation)
    set_maximal_velocities!(pbody,
        v=base_linear_velocity,
        ω=base_angular_velocity)

    previd = pbody.id
    for body in mechanism.bodies[2:end]
        set_maximal_configurations!(get_body(mechanism, previd), body,
            parent_vertex=vert12,
            child_vertex=vert11)
        set_maximal_velocities!(get_body(mechanism, previd), body,
            parent_vertex=vert12,
            child_vertex=vert11,
            Δv=relative_linear_velocity,
            Δω=relative_angular_velocity)
        previd = body.id
    end
end


text:
function flux(u, q, w, v, x, t, mu, eta)
    f = mu[1]*q + u[1]*mu[2:3];
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    x1 = x[1];
    x2 = x[2];
    v1 = mu[2];
    v2 = mu[3];
    s = (v1*x2*exp(v1)*exp(v2)*(exp(v2*(x2 - 1)) - 1)*(exp(v1*(x1 - 1)) + v1*x1*exp(v1*(x1 - 1)) - 1))/((exp(v1) - 1)*(exp(v2) - 1)) - (v2*x1*exp(v1 + v2*x2)*(exp(v1*(x1 - 1)) - 1)*(v2*x2 + 2))/((exp(v1) - 1)*(exp(v2) - 1)) - (v1*x2*exp(v2 + v1*x1)*(exp(v2*(x2 - 1)) - 1)*(v1*x1 + 2))/((exp(v1) - 1)*(exp(v2) - 1)) + (v2*x1*exp(v1)*exp(v2)*(exp(v1*(x1 - 1)) - 1)*(exp(v2*(x2 - 1)) + v2*x2*exp(v2*(x2 - 1)) - 1))/((exp(v1) - 1)*(exp(v2) - 1));
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = 0.0;
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);
    fb = f[1]*n[1] + f[2]*n[2] + tau[1]*(u[1]-uhat[1]);
    return fb;
end
function initu(x, mu, eta)
    u0 = 0.0;
    return u0;
end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
# LEARNINGS:
# - fields of a struct are compared via object equality (===)
ex = [
"light red bags contain 1 bright white bag, 2 muted yellow bags.",
"dark orange bags contain 3 bright white bags, 4 muted yellow bags.",
"bright white bags contain 1 shiny gold bag.",
"muted yellow bags contain 2 shiny gold bags, 9 faded blue bags.",
"shiny gold bags contain 1 dark olive bag, 2 vibrant plum bags.",
"dark olive bags contain 3 faded blue bags, 4 dotted black bags.",
"vibrant plum bags contain 5 faded blue bags, 6 dotted black bags.",
"faded blue bags contain no other bags.",
"dotted black bags contain no other bags.",
]

struct BagType
    Kind::String
    Contents::Dict
end

function parse_contents(s::String)
    if s[length(s)] == '.'
        s = chop(s)
    end

    out = Dict()
    if s == "no other bags"
        return out
    end

    parts = split(s, ",")
    for p in parts
        # p is something like "5 faded blue bags"
        words = split(p)
        cNum = parse(Int, words[1])
        cKind = words[2] * " " * words[3]
        out[cKind] = cNum
    end

    return out
end

function parse_rule(s::String)
    parts = split(s, " bags contain ")
    kind = parts[1]
    contents = parse_contents(String(parts[2]))

    out = BagType(kind, contents)
    return out
end

@assert parse_contents("5 faded blue bags, 6 dotted black bags.") == Dict(
    "faded blue" => 5,
    "dotted black" => 6,
)

@assert parse_contents("no other bags.") == Dict()

# parse rule
out = parse_rule("light red bags contain 1 bright white bag, 2 muted yellow bags.")
@assert out.Kind == "light red"
@assert out.Contents == Dict(
    "bright white" => 1,
    "muted yellow" => 2,
)

out = parse_rule("maroon tsunami bags contain no other bags.")
@assert out.Kind == "maroon tsunami"
@assert out.Contents == Dict()

function run(fname)
    # read file
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    bagTypes = map(parse_rule, lines)
    toSearch = []
    push!(toSearch, "shiny gold")
    searched = []
    result = 0
    while length(toSearch) > 0
        item = pop!(toSearch)
        push!(searched, item)

        # If it's in this bagType AND we haven't searched it yet
        for bag in bagTypes
            if haskey(bag.Contents, item)
                # check if parent bag's kind has been searched yet
                if ! in(bag.Kind, searched) && ! in(bag.Kind, toSearch)
                    push!(toSearch, bag.Kind)
                    result += 1
                end
            end
        end
    end

    println("Ran ", fname, " got ", result)
    return result
end

function count_bags_inside(bagMap, bag)
    contents = bagMap[bag].Contents
    if length(contents) == 0
        return 0
    else
        total = 0
        for k in keys(contents)
            val = contents[k]
            # add number of bags at this layer
            total += val
            r = count_bags_inside(bagMap, k)
            # add number of bags inside each of the bags at this layer
            total += (val * r)
        end
        return total
    end
end

function run2(fname)
    # read file
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    bagTypes = map(parse_rule, lines)
    bagMap = Dict()
    for b in bagTypes
        bagMap[b.Kind] = b
    end

    result = count_bags_inside(bagMap, "shiny gold")

    println("[part 2] Ran ", fname, " got ", result)
    return result
end

@assert run("7ex.txt") == 4
run("7.txt")

@assert run2("7ex2.txt") == 126
run2("7.txt")


text:
# YAML Tokens.
# Each token must include at minimum member "span::Span".
abstract type Token end

# The '%YAML' directive.
struct DirectiveToken <: Token
    span::Span
    name::AbstractString
    value::Union{Tuple, Nothing}
end

# '---'
struct DocumentStartToken <: Token
    span::Span
end

# '...'
struct DocumentEndToken <: Token
    span::Span
end

# The stream start
struct StreamStartToken <: Token
    span::Span
    encoding::AbstractString
end

# The stream end
struct StreamEndToken <: Token
    span::Span
end

#
struct BlockSequenceStartToken <: Token
    span::Span
end

#
struct BlockMappingStartToken <: Token
    span::Span
end

#
struct BlockEndToken <: Token
    span::Span
end

# '['
struct FlowSequenceStartToken <: Token
    span::Span
end

# '{'
struct FlowMappingStartToken <: Token
    span::Span
end

# ']'
struct FlowSequenceEndToken <: Token
    span::Span
end

# '}'
struct FlowMappingEndToken <: Token
    span::Span
end

# '?' or nothing (simple keys).
struct KeyToken <: Token
    span::Span
end

# ':'
struct ValueToken <: Token
    span::Span
end

# '-'
struct BlockEntryToken <: Token
    span::Span
end

# ','
struct FlowEntryToken <: Token
    span::Span
end

# '*anchor'
struct AliasToken <: Token
    span::Span
    value::AbstractString
end

# '&anchor'
struct AnchorToken <: Token
    span::Span
    value::AbstractString
end

# '!handle!suffix'
struct TagToken <: Token
    span::Span
    value
end

# A scalar.
struct ScalarToken <: Token
    span::Span
    value::AbstractString
    plain::Bool
    style::Union{Char, Nothing}
end


text:
function mass(u, q, w, v, x, t, mu, eta)
    m = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0];
    return m;
end
function flux(u, q, w, v, x, t, mu, eta)
    gam = mu[1];
    gam1 = gam - 1.0;
    r = u[1];
    ru = u[2];
    rv = u[3];
    rE = u[4];
    bx = u[5];
    by = u[6];
    phi = u[7];
    r1 = 1/r;
    uv = ru*r1;
    vv = rv*r1;
    E = rE*r1;
    q = 0.5*(uv*uv+vv*vv);
    b = (bx*bx+by*by)*0.5;
    p = gam1*(rE-r*q-b-0.5*phi^2);
    h = E+p*r1+b*r1;
    uvb = uv*bx+vv*by;
    f = [ru, ru*uv+p+b-bx*bx, rv*uv-by*bx, ru*h-uvb*bx+phi*bx, phi, -(vv*bx-by*uv), bx, rv, ru*vv-bx*by, rv*vv+p+b-by*by, rv*h-uvb*by+phi*by, -(uv*by-bx*vv), phi, by];
    f = reshape(f,(7,2));
    return f;
end
function source(u, q, w, v, x, t, mu, eta)
    s = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
    return s;
end
function ubou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    ub = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0];
    return ub;
end
function fbou(u, q, w, v, x, t, mu, eta, uhat, n, tau)
    f = flux(u, q, w, v, x, t, mu, eta);
    fb = f[:,1]*n[1] + f[:,2]*n[2] + tau[1]*(u-uhat);
    return fb;
end
function initu(x, mu, eta)
    gam = mu[1];
    x1 = x[1];
    x2 = x[2];
    r = sqrt(x1*x1+x2*x2);
    k = 1/(2*pi);
    u01 = 1;
    u02 = u01-k*x2*exp(0.5*(1-r^2));
    u03 = u01+k*x1*exp(0.5*(1-r^2));
    u05 = -k*x2*exp(0.5*(1-r^2));
    u06 = k*x1*exp(0.5*(1-r^2));
    p = 1 + 1/(4*0.5)*(k^2*(1-2*0.5*r^2)-k^2*u01)*exp(2*0.5*(1-r^2));
    b = u06*u06 + u05*u05;
    v = u02*u02 + u03*u03;
    u04 = p/(gam-1) + u01/2*v + 0.5*b;
    u07 = 0;
    u0 = [u01, u02, u03, u04, u05, u06, u07];
    return u0;
end


text:
#///////////////////////////////////////
#// File Name: nelder_mead_bilevel_optimization_test.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2020/11/05
#// Description: Test code for src/nelder_mead_bilevel_optimization.jl
#///////////////////////////////////////

using Test
using LinearAlgebra

@testset "Nelder-Mead Simplex Bilevel Optimization Test" begin
    f(x, u) = x.^1.3 + u.^1.5;
    c(k, x, u) = sum(x.^2.5 + u.^2.5);
    h(x) = 1.0;
    W(k) = Matrix(0.01I, 2, 2);
    N = 10;

    x_0 = zeros(2)
    u_array = [0.1*ones(2) for ii = 1 : N];

    problem = FiniteHorizonRiskSensitiveOptimalControlProblem(f, c, h, W, N);
    kl_bound = 1.0;
    nm_solver = NelderMeadBilevelOptimizationSolver(iter_max=20, ϵ=1e-3, θ_high_init=10.0, θ_low_init=1e-8);

    θ_opt, x_array, l_array, L_array, c_opt =
    solve!(nm_solver, problem, x_0, u_array, kl_bound=kl_bound, verbose=false)
    @test !isinf(c_opt)
    @test !isnan(θ_opt)
    c_low_init = compute_cost_worker(nm_solver, problem, x_0, u_array, nm_solver.θ_low_init, kl_bound)
    c_high_init = compute_cost_worker(nm_solver, problem, x_0, u_array, nm_solver.θ_high_init, kl_bound)
    @test !isinf(c_low_init) && !isinf(c_high_init);
    @test c_opt <= c_low_init && c_opt <= c_high_init
end


text:
using BenchmarkTools
using Random
using ReactiveMP

import Pkg;

# Pkg.resolve()
Pkg.instantiate()
Pkg.status()

# Define a parent BenchmarkGroup to contain our suite
const SUITE = BenchmarkGroup()

SUITE["models"] = BenchmarkGroup([ "models", "ssm", "graphppl" ])

# Simple Linear Gaussian State Space Model Benchmarks
# ------------------------------------------------------------ #
include("models/lgssm1.jl")

SUITE["models"]["lgssm1"] = BenchmarkGroup([ "linear", "gaussian", "ssm", "univariate" ])

for lgssm1_size in [ 100, 500 ]
    # Model creation benchmark
    SUITE["models"]["lgssm1"]["creation_$lgssm1_size"]  = @benchmarkable LGSSM1Benchmark.lgssm($lgssm1_size)
    # Inference benchmark
    SUITE["models"]["lgssm1"]["inference_$lgssm1_size"] = @benchmarkable LGSSM1Benchmark.benchmark(input) setup=(input=LGSSM1Benchmark.generate_input(MersenneTwister(1234), $lgssm1_size))
end
# ------------------------------------------------------------ #

# Multivariate Linear Gaussian State Space Model Benchmarks
# ------------------------------------------------------------ #
include("models/lgssm2.jl")

SUITE["models"]["lgssm2"] = BenchmarkGroup([ "linear", "gaussian", "ssm", "multivariate" ])

for lgssm2_size in [ 100, 500 ]
    # Model creation benchmark
    SUITE["models"]["lgssm2"]["creation_$lgssm2_size"]  = @benchmarkable LGSSM2Benchmark.lgssm($lgssm2_size)
    # Inference benchmark
    SUITE["models"]["lgssm2"]["inference_$lgssm2_size"] = @benchmarkable LGSSM2Benchmark.benchmark(input) setup=(input=LGSSM2Benchmark.generate_input(MersenneTwister(1234), $lgssm2_size))
end
# ------------------------------------------------------------ #

# Hierarchical Gaussian Filter Model Benchmarks
# ------------------------------------------------------------ #
include("models/hgf1.jl")

SUITE["models"]["hgf1"] = BenchmarkGroup([ "hierarchical", "gaussian", "ssm", "univariate" ])

SUITE["models"]["hgf1"]["creation"]  = @benchmarkable HGF1Benchmark.hgf()

for hgf1_size in [ 100, 500, 1000 ]
    # Inference benchmark
    SUITE["models"]["hgf1"]["inference_$hgf1_size"] = @benchmarkable HGF1Benchmark.benchmark(input) setup=(input=HGF1Benchmark.generate_input(MersenneTwister(1234), $hgf1_size))
end
# ------------------------------------------------------------ #

# Hidden Markov Model Benchmarks
# ------------------------------------------------------------ #
include("models/hmm1.jl")

SUITE["models"]["hmm1"] = BenchmarkGroup([ "hmm", "ssm", "discrete" ])

for hmm1_size in [ 100, 500 ]
    # Creation Benchmark
    SUITE["models"]["hmm1"]["creation_$hmm1_size"]  = @benchmarkable HMM1Benchmark.hmm($hmm1_size)
    # Inference benchmark
    SUITE["models"]["hmm1"]["inference_$hmm1_size"] = @benchmarkable HMM1Benchmark.benchmark(input) setup=(input=HMM1Benchmark.generate_input(MersenneTwister(1234), $hmm1_size))
end
# ------------------------------------------------------------ #

BenchmarkTools.warmup(SUITE)


text:
using JLD2
using DSGE
m = KrusellSmithCT()
T, R, ~, inverse_basis = solvect(m)

n_vars = get_setting(m, :n_vars)
n_expectation_errors = get_setting(m, :n_expectation_errors)
n_shocks = get_setting(m, :n_shocks)
x = zeros(Float64, n_vars)

nst = size(T,1) ###
C = zeros(nst)

n = 40 #100
s = zeros(nst,n) #s = zeros(n,406)
s[:,1] = inverse_basis' * x
e = randn(n)

dt = 1/2

for i = 2:n
    s[:, i] = s[:, i-1] + T*s[:,i-1]*dt + R*e[i]*(sqrt(dt))
end
data = zeros(1,Int(n/2))
data[1, :] = [s[nst - 2, 2 * i] for i = 1:Int(n/2)]

save("simulated_data.jld2", "data", data)


text:
function default_attributes(::Type{LAxis})
    Attributes(
        xlabel = "x label",
        ylabel = "y label",
        title = "Title",
        titlefont = "DejaVu Sans",
        titlesize = 30f0,
        titlegap = 10f0,
        titlevisible = true,
        titlealign = :center,
        xlabelcolor = RGBf0(0, 0, 0),
        ylabelcolor = RGBf0(0, 0, 0),
        xlabelsize = 20f0,
        ylabelsize = 20f0,
        xlabelvisible = true,
        ylabelvisible = true,
        xlabelpadding = 5f0,
        ylabelpadding = 5f0,
        xticklabelsize = 20f0,
        yticklabelsize = 20f0,
        xticklabelsvisible = true,
        yticklabelsvisible = true,
        xticklabelspace = 20f0,
        yticklabelspace = 50f0,
        xticklabelpad = 5f0,
        yticklabelpad = 5f0,
        xticklabelrotation = 0f0,
        yticklabelrotation = 0f0,
        xticklabelalign = (:center, :top),
        yticklabelalign = (:right, :center),
        xticksize = 10f0,
        yticksize = 10f0,
        xticksvisible = true,
        yticksvisible = true,
        xtickalign = 0f0,
        ytickalign = 0f0,
        xtickwidth = 1f0,
        ytickwidth = 1f0,
        xtickcolor = RGBf0(0, 0, 0),
        ytickcolor = RGBf0(0, 0, 0),
        xpanlock = false,
        ypanlock = false,
        xzoomlock = false,
        yzoomlock = false,
        spinewidth = 1f0,
        xgridvisible = true,
        ygridvisible = true,
        xgridwidth = 1f0,
        ygridwidth = 1f0,
        xgridcolor = RGBAf0(0, 0, 0, 0.1),
        ygridcolor = RGBAf0(0, 0, 0, 0.1),
        xgridstyle = nothing,
        ygridstyle = nothing,
        xspinevisible = true,
        yspinevisible = true,
        xoppositespinevisible = true,
        yoppositespinevisible = true,
        xspinecolor = :black,
        yspinecolor = :black,
        xoppositespinecolor = :black,
        yoppositespinecolor = :black,
        aspect = nothing,
        valign = :center,
        halign = :center,
        width = nothing,
        height = nothing,
        maxsize = (Inf32, Inf32),
        xautolimitmargin = (0.05f0, 0.05f0),
        yautolimitmargin = (0.05f0, 0.05f0),
        xticks = AutoLinearTicks(100f0),
        yticks = AutoLinearTicks(100f0),
        panbutton = AbstractPlotting.Mouse.right,
        xpankey = AbstractPlotting.Keyboard.x,
        ypankey = AbstractPlotting.Keyboard.y,
        xzoomkey = AbstractPlotting.Keyboard.x,
        yzoomkey = AbstractPlotting.Keyboard.y,
        xaxisposition = :bottom,
        yaxisposition = :left,
        xtrimspine = false,
        ytrimspine = false,
    )
end

function default_attributes(::Type{LColorbar})
    Attributes(
        label = "label",
        labelcolor = RGBf0(0, 0, 0),
        labelsize = 20f0,
        labelvisible = true,
        labelpadding = 5f0,
        ticklabelsize = 20f0,
        ticklabelsvisible = true,
        ticksize = 10f0,
        ticksvisible = true,
        ticklabelspace = 30f0,
        ticklabelpad = 5f0,
        tickalign = 0f0,
        tickwidth = 1f0,
        tickcolor = RGBf0(0, 0, 0),
        ticklabelalign = (:left, :center),
        spinewidth = 1f0,
        idealtickdistance = 100f0,
        topspinevisible = true,
        rightspinevisible = true,
        leftspinevisible = true,
        bottomspinevisible = true,
        topspinecolor = RGBf0(0, 0, 0),
        leftspinecolor = RGBf0(0, 0, 0),
        rightspinecolor = RGBf0(0, 0, 0),
        bottomspinecolor = RGBf0(0, 0, 0),
        valign = :center,
        halign = :center,
        vertical = true,
        flipaxisposition = true,
        width = nothing,
        height = nothing,
        colormap = :viridis,
        limits = (0f0, 1f0),
    )
end

function default_attributes(::Type{LText})
    Attributes(
        text = "Text",
        visible = true,
        color = RGBf0(0, 0, 0),
        textsize = 20f0,
        font = "Dejavu Sans",
        valign = :center,
        halign = :center,
        rotation = 0f0,
        padding = (0f0, 0f0, 0f0, 0f0),
        height = Auto(),
        width = Auto(),
    )
end

function default_attributes(::Type{LRect})
    Attributes(
        visible = true,
        color = RGBf0(0.9, 0.9, 0.9),
        valign = :center,
        halign = :center,
        padding = (0f0, 0f0, 0f0, 0f0),
        strokewidth = 2f0,
        strokevisible = true,
        strokecolor = RGBf0(0, 0, 0),
        width = nothing,
        height = nothing,
    )
end

function default_attributes(::Type{LButton})
    Attributes(
        halign = :center,
        valign = :center,
        padding = (10f0, 10f0, 10f0, 10f0),
        textsize = 20f0,
        label = "Button",
        font = "Dejavu Sans",
        width = Auto(true),
        height = Auto(true),
        cornerradius = 4,
        cornersegments = 10,
        strokewidth = 2f0,
        strokecolor = :transparent,
        buttoncolor = RGBf0(0.9, 0.9, 0.9),
        labelcolor = :black,
        labelcolor_hover = :black,
        labelcolor_active = :white,
        buttoncolor_active = COLOR_ACCENT[],
        # buttoncolor_hover = RGBf0(0.8, 0.8, 0.8),
        buttoncolor_hover = COLOR_ACCENT_DIMMED[],
        clicks = 0,
    )
end

function default_attributes(::Type{AxisContent})
    Attributes(
        xautolimit = true,
        yautolimit = true,
    )
end

function default_attributes(::Type{LineAxis})
    Attributes(
        endpoints = (Point2f0(0, 0), Point2f0(100, 0)),
        trimspine = false,
        limits = (0f0, 100f0),
        flipped = false,
        ticksize = 10f0,
        tickwidth = 1f0,
        tickcolor = RGBf0(0, 0, 0),
        tickalign = 0f0,
        ticks = AutoLinearTicks(100f0),
        ticklabelalign = (:center, :top),
        ticksvisible = true,
        ticklabelrotation = 0f0,
        ticklabelsize = 20f0,
        ticklabelsvisible = true,
        spinewidth = 1f0,
        label = "label",
        labelsize = 20f0,
        labelcolor = RGBf0(0, 0, 0),
        labelvisible = true,
        ticklabelspace = 30f0,
        ticklabelpad = 5f0,
        labelpadding = 10f0,
    )
end

function default_attributes(::Type{LSlider})
    Attributes(
        linewidth = 4f0,
        halign = :center,
        valign = :center,
        # vertical = true,
        width = nothing,
        height = Auto(true),
        range = 0:10,
        buttonradius = 7f0,
        startvalue = 0,
        value = 0,
        color_active_dimmed = COLOR_ACCENT_DIMMED[],
        color_active = COLOR_ACCENT[],
        color_inactive = RGBf0(0.9, 0.9, 0.9),
        buttoncolor_inactive = RGBf0(1, 1, 1),
        horizontal = true,
        buttonstrokewidth = 4f0,
    )
end

function default_attributes(::Type{LToggle})
    Attributes(
        halign = :center,
        valign = :center,
        width = 60,
        height = 30,
        cornersegments = 10,
        # strokewidth = 2f0,
        # strokecolor = :transparent,
        framecolor_inactive = RGBf0(0.9, 0.9, 0.9),
        framecolor_active = COLOR_ACCENT[],
        # buttoncolor = RGBf0(0.2, 0.2, 0.2),
        buttoncolor = RGBf0(1, 1, 1),
        active = false,
        toggleduration = 0.2,
        rimfraction = 0.25,
    )
end

function default_attributes(::Type{LLegend})
    Attributes(
        halign = :center,
        valign = :center,
        width = Auto(true),
        height = Auto(false),
        title = "Legend",
        titlefont = "Dejavu Sans",
        titlesize = 20f0,
        titlealign = :center,
        titlevisible = true,
        labelsize = 20f0,
        labelfont = "Dejavu Sans",
        labelcolor = :black,
        labelhalign = :left,
        labelvalign = :center,
        padding = (10f0, 10f0, 10f0, 10f0),
        margin = (0f0, 0f0, 0f0, 0f0),
        bgcolor = :white,
        strokecolor = :black,
        strokewidth = 1f0,
        patchsize = (40f0, 40f0),
        patchstrokecolor = :transparent,
        patchstrokewidth = 1f0,
        patchcolor = RGBf0(0.97, 0.97, 0.97),
        label = "undefined",
        ncols = 1,
        colgap = 20,
        rowgap = 4,
        patchlabelgap = 5,
        linepoints = [Point2f0(0, 0.5), Point2f0(1, 0.5)],
        linewidth = 3,
        markerpoints = [Point2f0(0.5, 0.5)],
        markersize = 16,
        markerstrokewidth = 2,
        polypoints = [Point2f0(0.2, 0.2), Point2f0(0.8, 0.2), Point2f0(0.8, 0.8), Point2f0(0.2, 0.8)],
        polystrokewidth = 2,
    )
end

function attributenames(::Type{LegendEntry})
    (:label, :labelsize, :labelfont, :labelcolor, :labelhalign, :labelvalign,
        :patchsize, :patchstrokecolor, :patchstrokewidth, :patchcolor,
        :linepoints, :markerpoints, :markersize, :markerstrokewidth, :linewidth,
        :polypoints, :polystrokewidth)
end

function extractattributes(attributes::Attributes, typ::Type)
    extracted = Attributes()
    for name in attributenames(typ)
        if haskey(attributes, name)
            extracted[name] = attributes[name]
        end
    end
    extracted
end

function default_attributes(::Type{GridLayout})
    Attributes(
        halign = :center,
        valign = :center,
        width = Auto(),
        height = Auto(),
    )
end


text:
## Define the new measure evaluation method
# Make alias for our new method
const NewEvalMethod = :NewEvalMethod # REPLACE WITH ACTUAL

# Extend generate_supports_and_coeffs for scalar params
function InfiniteOpt.MeasureEvalMethods.generate_supports_and_coeffs(
    set::InfiniteOpt.IntervalSet,
    params::Union{InfiniteOpt.ParameterRef, AbstractArray{<:InfiniteOpt.ParameterRef}},
    num_supports::Int,
    lb::Number,
    ub::Number,
    method::Val{NewEvalMethod})::Tuple # REPLACE WITH ACTUAL ALIAS
    # REPLACE WITH ACTUAL FUNCTIONALITY
    increment = (ub - lb) / (num_supports - 1)
    supports = [lb + (i - 1) * increment for i in 1:num_supports]
    # MAKE SURE IT RETURNS APPROPRIATE DATA IN THE TUPLE
    return (supports, ones(num_supports) / num_supports * (ub - lb))
end

# Extend generate_supports_and_coeffs for multi-dimensional params
function InfiniteOpt.MeasureEvalMethods.generate_supports_and_coeffs(
    set::InfiniteOpt.IntervalSet,
    params::Union{InfiniteOpt.ParameterRef, AbstractArray{<:InfiniteOpt.ParameterRef}},
    num_supports::Int,
    lb::JuMPC.SparseAxisArray,
    ub::JuMPC.SparseAxisArray,
    method::Val{NewEvalMethod}; independent::Bool = true)::Tuple # REPLACE WITH ACTUAL ARGUMENTS
    # REPLACE WITH ACTUAL FUNCTIONALITY
    if !independent
        @warn("The method is implemented for independent multivariate parameters.")
    end
    supports_dict = Dict()
    for i in eachindex(lb)
        (supports_dict[i], _) = generate_supports_and_coeffs(set, params,
                                 num_supports, lb[i], ub[i], Val(NewEvalMethod))
    end
    supports = Array{JuMPC.SparseAxisArray, 1}(undef, num_supports)
    for j in 1:num_supports
        supports[j] = JuMP.Containers.SparseAxisArray(Dict(k => supports_dict[k][j]
                                                        for k in eachindex(lb)))
    end
    # MAKE SURE IT RETURNS APPROPRIATE DATA IN THE TUPLE
    return (supports, ones(num_supports) / num_supports * prod(ub .- lb))
end


text:
function remove_fully_collected!(id::Int, set::Symbol, rarity::Int; openable_cards)
    filter!(z -> z != id, openable_cards.by_rarity[rarity])
    if haskey(openable_cards.by_set_rarity, (set,rarity))
        filter!(z -> z != id, openable_cards.by_set_rarity[(set,rarity)])
    end

    openable_cards
end

function remove_fully_collected!(id::Int, rarity::Int; openable_cards)
    for set in openable_cards.sets
        remove_fully_collected!(id, set, rarity; openable_cards = openable_cards)
    end

    openable_cards
end

function remove_fully_collected!(id::Int; openable_cards)
    for r in 1:4
        remove_fully_collected!(id, r; openable_cards = openable_cards)
    end

    openable_cards
end

max_stars_repo_name: asieev/ArenaSim

max_issues_repo_name: asieev/ArenaSim

max_forks_repo_name: asieev/ArenaSim


text:
include("setup.jl")

NON_STDLIB_TESTS &&
@testset "ChaosBufferStream" begin
    @testset "constant usage" begin
        io = BufferStream()
        cio = ChaosBufferStream(io; chunksizes=[17], sleepamnts=[0.001])
        write(io, rand(UInt8, 30))
        close(io)

        # Test that data comes out in 17-byte chunks (except for the last)
        buff = Array{UInt8}(undef, 30)
        t = @elapsed begin
            @test readbytes!(cio, buff, 30) == 17
            @test readbytes!(cio, buff, 30) == 13
        end
        @test t >= 0.001
    end

    @testset "random usage" begin
        io = BufferStream()
        chunksizes = 5:10
        cio = ChaosBufferStream(io; chunksizes=chunksizes, sleepamnts=[0.0])
        write(io, rand(UInt8, 3000))
        close(io)

        buff = Array{UInt8}(undef, 10)
        while !eof(cio)
            r = readbytes!(cio, buff, 10)
            # In normal operation, the chunk size must be one of
            # the given chunksizes, but at the end of the stream
            # it is allowed to be less.
            if !eof(cio)
                @test r ∈ chunksizes
            else
                @test r <= maximum(chunksizes)
            end
        end
    end
end

@testset "empty tarball" begin
    @testset "empty file as tarball" begin
        tarball = devnull
        @test Tar.list(tarball) == Tar.Header[]
        skel = tempname()
        dir = Tar.extract(tarball, skeleton=skel)
        @test isempty(readdir(dir))
        rm(dir, recursive=true)
        @test isfile(skel)
        @test Tar.list(skel) == Tar.Header[]
        @test Tar.list(skel, raw=true) == Tar.Header[]
        rm(skel)
    end

    @testset "create an empty tarball" begin
        dir = mktempdir()
        tarball = Tar.create(dir)
        rm(dir, recursive=true)
        @test Tar.list(tarball) == [Tar.Header(".", :directory, 0o755, 0, "")]
        @test Tar.list(tarball, raw=true) == [Tar.Header(".", :directory, 0o755, 0, "")]
        test_empty_hashes(tarball)
        skel = tempname()
        dir = Tar.extract(tarball, skeleton=skel)
        @test isempty(readdir(dir))
        rm(dir, recursive=true)
        @test isfile(skel)
        @test Tar.list(skel) == [Tar.Header(".", :directory, 0o755, 0, "")]
        @test Tar.list(skel, raw=true) == [Tar.Header(".", :directory, 0o755, 0, "")]
        rm(skel)
        open(tarball, append=true) do io
            write(io, zeros(UInt8, 512))
        end
        test_empty_hashes(tarball)
        dir = Tar.extract(tarball)
        @test isempty(readdir(dir))
        rm(dir, recursive=true)
    end
end

@testset "test tarball" begin
    tarball, hash = make_test_tarball()
    @testset "Tar.tree_hash" begin
        arg_readers(tarball) do tar
            @arg_test tar @test Tar.tree_hash(tar, skip_empty=true) == hash
            @arg_test tar @test empty_tree_sha1 == Tar.tree_hash(hdr->false, tar)
            @arg_test tar @test empty_tree_sha1 ==
                Tar.tree_hash(hdr->false, tar, algorithm="git-sha1")
            @arg_test tar @test empty_tree_sha256 ==
                Tar.tree_hash(hdr->false, tar, algorithm="git-sha256")
        end
    end
    @testset "Tar.list & check properties" begin
        headers = Tar.list(tarball)
        @test issorted(headers, by = hdr -> hdr.path)
        for hdr in headers
            @test !isempty(hdr.path)
            @test hdr.path[end] != '/'
            @test hdr.type in (:file, :directory, :symlink)
            if hdr.type == :file
                @test hdr.mode in (0o644, 0o755)
                @test isempty(hdr.link)
            elseif hdr.type == :directory
                @test hdr.mode == 0o755
                @test hdr.size == 0
                @test isempty(hdr.link)
            elseif hdr.type == :symlink
                @test hdr.mode == 0o755
                @test hdr.size == 0
                @test !isempty(hdr.link)
            end
        end
        @testset "Tar.list from IO, process, pipeline" begin
            arg_readers(tarball) do tar
                @arg_test tar begin
                    @test headers == Tar.list(tar)
                end
            end
        end
    end
    if @isdefined(gtar)
        @testset "extract with `tar` command" begin
            root = mktempdir()
            gtar(gtar -> run(`$gtar -C $root -xf $tarball`))
            check_tree_hash(hash, root)
        end
    end
    @testset "Tar.extract" begin
        arg_readers(tarball) do tar
            @arg_test tar begin
                root = Tar.extract(tar)
                check_tree_hash(hash, root)
            end
        end
    end
    open(tarball, append=true) do io
        write(io, zeros(UInt8, 512))
    end
    @testset "Tar.extract with trailing zeros" begin
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
    rm(tarball)
end

@testset "truncated tarball" begin
    tarball, hash = make_test_tarball()
    open(tarball, "a") do io
        truncate(io, div(filesize(tarball),2))
    end

    @testset "tree_hash" begin
        # Ensure that this throws
        @test_throws EOFError Tar.tree_hash(tarball)
    end
end

if @isdefined(gtar)
    @testset "POSIX extended headers" begin
        # make a test POSIX tarball with GNU `tar` from Tar_jll instead of Tar.create
        tarball, hash = make_test_tarball() do root
            tarball, io = mktemp()
            Tar.write_extended_header(io, type = :g, ["comment" => "Julia Rocks!"])
            close(io)
            gtar(gtar -> run(`$gtar --format=posix -C $root --append -f $tarball .`))
            return tarball
        end
        # TODO: check that extended headers contain `mtime` etc.
        @test Tar.tree_hash(tarball, skip_empty=true) == hash
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
    @testset "GNU extensions" begin
        # make a test GNU tarball with GNU `tar` from Tar_jll instead of Tar.create
        tarball, hash = make_test_tarball() do root
            tarball = tempname()
            gtar(gtar -> run(`$gtar --format=gnu -C $root -cf $tarball .`))
            return tarball
        end
        hdrs = Tar.list(tarball, raw=true)
        # test that both long link and long name entries are created
        @test any(h.path == "././@LongLink" && h.type == :L for h in hdrs)
        @test any(h.path == "././@LongLink" && h.type == :K for h in hdrs)
        # test that Tar can extract these GNU entries correctly
        @test Tar.tree_hash(tarball, skip_empty=true) == hash
        root = Tar.extract(tarball)
        check_tree_hash(hash, root)
    end
end

@testset "symlink attacks" begin
    # not dangerous but still not allowed
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("dir", :directory, 0o755, 0, ""))
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "dir"))
    Tar.write_header(io, Tar.Header("link/target", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    @test_throws ErrorException Tar.tree_hash(tarball)
    rm(tarball)
    # attempt to write through relative link out of root
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, "../target"))
    Tar.write_header(io, Tar.Header("link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    @test_throws ErrorException Tar.tree_hash(tarball)
    rm(tarball)
    # attempt to write through absolute link
    tmp = mktempdir()
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, tmp))
    Tar.write_header(io, Tar.Header("link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    @test_throws ErrorException Tar.tree_hash(tarball)
    rm(tarball)
    # same attack with some obfuscation
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, tmp))
    Tar.write_header(io, Tar.Header("./link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    @test_throws ErrorException Tar.tree_hash(tarball)
    rm(tarball)
    # same attack with different obfuscation
    tarball, io = mktemp()
    Tar.write_header(io, Tar.Header("link", :symlink, 0o755, 0, tmp))
    Tar.write_header(io, Tar.Header("dir/../link/attack", :file, 0o644, 0, ""))
    close(io)
    @test_throws ErrorException Tar.extract(tarball)
    @test_throws ErrorException Tar.tree_hash(tarball)
    rm(tarball)
    # check temp dir is empty, remove it
    @test isempty(readdir(tmp))
    rm(tmp)
end

!Sys.iswindows() &&
@testset "symlink overwrite" begin
    tmp = mktempdir()
    @testset "allow overwriting a symlink" begin
        tarball₁, io = mktemp()
        Tar.write_header(io, Tar.Header("path", :symlink, 0o755, 0, tmp))
        Tar.write_header(io, Tar.Header("path", :file, 0o644, 0, ""))
        close(io)
        hash = Tar.tree_hash(tarball₁)
        tree₁ = Tar.extract(tarball₁)
        @test hash == tree_hash(tree₁)
        tarball₂, io = mktemp()
        Tar.write_header(io, Tar.Header("path", :file, 0o644, 0, ""))
        close(io)
        @test hash == Tar.tree_hash(tarball₂)
        tree₂ = Tar.extract(tarball₂)
        @test hash == tree_hash(tree₂)
        rm(tree₁, recursive=true)
        rm(tree₂, recursive=true)
        rm(tarball₁)
        rm(tarball₂)
    end
    @testset "allow write into directory overwriting a symlink" begin
        # make sure "path" is removed from links set
        tarball₁, io = mktemp()
        Tar.write_header(io, Tar.Header("path", :symlink, 0o755, 0, tmp))
        Tar.write_header(io, Tar.Header("path", :directory, 0o755, 0, ""))
        Tar.write_header(io, Tar.Header("path/file", :file, 0o644, 0, ""))
        close(io)
        hash = Tar.tree_hash(tarball₁)
        tree₁ = Tar.extract(tarball₁)
        @test hash == tree_hash(tree₁)
        tarball₂, io = mktemp()
        Tar.write_header(io, Tar.Header("path/file", :file, 0o644, 0, ""))
        close(io)
        @test hash == Tar.tree_hash(tarball₂)
        tree₂ = Tar.extract(tarball₂)
        @test hash == tree_hash(tree₂)
        rm(tree₁, recursive=true)
        rm(tree₂, recursive=true)
        rm(tarball₁)
        rm(tarball₂)
    end
    # check temp dir is empty, remove it
    @test isempty(readdir(tmp))
    rm(tmp)
end

@testset "copy symlinks" begin
    tmp = mktempdir()
    data₁ = randstring(12)
    data₂ = randstring(12)
    tarball, io = mktemp()
    tar_write_file(io, "file", data₁)
    tar_write_link(io, "link-file", "file")
    tar_write_link(io, "link-file-slash", "file/")
    tar_write_link(io, "link-file-slash-dot", "file/.")
    tar_write_link(io, "link-dot-slash-file", "./file")
    tar_write_link(io, "link-dir-dot-dot-file", "dir/../file")
    tar_write_link(io, "link-non-dot-dot-file", "non/../file")
    tar_write_link(io, "link-dir-non-dot-dot-dot-dot-file", "dir/non/../../file")
    tar_write_link(io, "link-self", "link-self")
    tar_write_link(io, "link-cycle-A", "link-cycle-B")
    tar_write_link(io, "link-cycle-B", "link-cycle-A")
    tar_write_link(io, "link-cycle-C", "link-cycle-B")
    tar_write_dir(io,  "dir")
    tar_write_link(io, "link-tmp", tmp)
    tar_write_link(io, "link-dot-dot", "..")
    tar_write_link(io, "link-dot", ".")
    tar_write_link(io, "link-dir", "dir")
    tar_write_link(io, "link-dir-slash", "dir/")
    tar_write_link(io, "link-dir-slash-dot", "dir/.")
    tar_write_link(io, "link-dot-slash-dir", "./dir")
    tar_write_link(io, "link-dot-slash-dir-slash", "./dir/")
    tar_write_link(io, "link-dot-slash-dir-slash-dot", "./dir/.")
    tar_write_link(io, "dir/link-file", "file")
    tar_write_link(io, "dir/link-dot-dot-dir-file", "../dir/file")
    tar_write_link(io, "dir/link-dot-dot-link-dir-file", "../link-dir/file")
    tar_write_link(io, "dir/link-dot-dot-file", "../file")
    tar_write_link(io, "dir/link-dot", ".")
    tar_write_link(io, "dir/link-dot-dot", "..")
    tar_write_link(io, "dir/link-dot-dot-dir", "../dir")
    tar_write_link(io, "dir/link-dot-dot-dir-self", "../dir/link-self")
    tar_write_file(io, "dir/file", data₂)
    close(io)
    # some test utilities (capturing `dir` defined below)
    test_none(path::String) = @test !ispath(joinpath(dir, path))
    test_dir(path::String) = @test isdir(joinpath(dir, path))
    function test_dir(a::String, b::String)
        A = joinpath(dir, a)
        B = joinpath(dir, b)
        @test isdir(A)
        @test isdir(B)
        @test read(Tar.create(A)) == read(Tar.create(B))
    end
    function test_file(path::String, data::String)
        path = joinpath(dir, path)
        @test isfile(path)
        @test read(path, String) == data
    end
    dir = Tar.extract(tarball, copy_symlinks=true)
    test_file("file", data₁)
    test_file("link-file", data₁)
    test_none("link-file-slash")
    test_none("link-file-slash-dot")
    test_file("link-dot-slash-file", data₁)
    test_file("link-dir-dot-dot-file", data₁)
    test_none("link-non-dot-dot-file")
    test_none("link-dir-non-dot-dot-dot-dot-file")
    test_none("link-cycle-A")
    test_none("link-cycle-B")
    test_none("link-cycle-C")
    test_dir("dir")
    test_none("link-tmp")
    test_none("link-dot-dot")
    test_none("link-dot")
    test_dir("link-dir", "dir")
    test_dir("link-dir-slash", "dir")
    test_dir("link-dir-slash-dot", "dir")
    test_dir("link-dot-slash-dir", "dir")
    test_dir("link-dot-slash-dir-slash", "dir")
    test_dir("link-dot-slash-dir-slash-dot", "dir")
    test_file("dir/link-file", data₂)
    test_file("dir/link-dot-dot-dir-file", data₂)
    test_file("dir/link-dot-dot-link-dir-file", data₂)
    test_file("dir/link-dot-dot-file", data₁)
    test_none("dir/link-dot")
    test_none("dir/link-dot-dot")
    test_none("dir/link-dot-dot-dir")
    test_none("dir/link-dot-dot-dir-self")
    test_file("dir/file", data₂)
    rm(dir, recursive=true)
    # check temp dir is empty, remove it
    @test isempty(readdir(tmp))
    rm(tmp)
end

@testset "API: create" begin
    local bytes

    @testset "without predicate" begin
        dir = make_test_dir()
        @test !any(splitext(name)[2] == ".skip" for name in readdir(dir))

        # create(dir)
        tarball = Tar.create(dir)
        bytes = read(tarball)
        @test isfile(tarball)
        rm(tarball)

        # create(dir, tarball)
        arg_writers() do tarball, tar
            @arg_test tar begin
                @test tar == Tar.create(dir, tar)
            end
            @test read(tarball) == bytes
        end

        # cleanup
        rm(dir, recursive=true)
    end

    @testset "with predicate" begin
        dir = make_test_dir(true)
        @test any(splitext(name)[2] == ".skip" for name in readdir(dir))
        predicate = path -> splitext(path)[2] != ".skip"

        # create(predicate, dir)
        tarball = Tar.create(predicate, dir)
        @test read(tarball) == bytes
        rm(tarball)

        # create(predicate, dir, tarball)
        arg_writers() do tarball, tar
            @arg_test tar begin
                @test tar == Tar.create(predicate, dir, tar)
            end
            @test read(tarball) == bytes
        end

        # cleanup
        rm(dir, recursive=true)
    end

    # In this issue we've seen that symlinking a directory caused files inside
    # the directory to become read-only.  Guard against Tar.jl doing something
    # weird like that.
    @testset "Issue Pkg#2185" begin
        mktempdir() do dir
            root = joinpath(dir, "root")
            target = joinpath("lib", "icu", "67.1")
            link = joinpath("lib", "icu", "current")
            file = joinpath(target, "file")
            dir_mode = 0o755
            file_mode = 0o644
            mkpath(joinpath(root, target))
            touch(joinpath(root, file))
            chmod(joinpath(root, file), dir_mode)
            chmod(joinpath(root, file), file_mode)
            symlink(basename(target), joinpath(root, link))
            tarball = Tar.create(root, joinpath(dir, "test.tar"))
            files = Tar.list(tarball)
            # Make sure the file and the symlink have the expected permissions.
            # Note: in old versions of Julia, the file has always permission 755 on Windows
            @test Tar.Header(replace(file, "\\" => "/"), :file, VERSION ≤ v"1.6.0-DEV.1683" && Sys.iswindows() ? 0o755 : file_mode, 0, "") in files
            @test Tar.Header(replace(link, "\\" => "/"), :symlink, dir_mode, 0, basename(target)) in files
        end
    end
end

@testset "API: list" begin
    dir = make_test_dir()
    tarball = Tar.create(dir)
    rm(dir, recursive=true)
    n = length(test_dir_paths)

    # list([predicate,] tarball)
    arg_readers(tarball) do tar
        @arg_test tar begin
            headers = Tar.list(tar)
            @test test_dir_paths == [hdr.path for hdr in headers]
        end
        @arg_test tar @test n == tar_count(tar)
    end

    # add a sketchy entry to tarball
    open(tarball, append=true) do io
        Tar.write_header(io, Tar.Header("/bad", :file, 0o644, 0, ""))
    end
    paths = push!(copy(test_dir_paths), "/bad")

    # list([predicate,] tarball; strict=true|false)
    arg_readers(tarball) do tar
        @arg_test tar @test_throws ErrorException Tar.list(tar)
        @arg_test tar @test_throws ErrorException Tar.list(tar, strict=true)
        @arg_test tar begin
            headers = Tar.list(tar, strict=false)
            @test paths == [hdr.path for hdr in headers]
        end
        @arg_test tar @test_throws ErrorException tar_count(tar)
        @arg_test tar @test_throws ErrorException tar_count(tar, strict=true)
        @arg_test tar @test n + 1 == tar_count(tar, strict=false)
    end
end

@testset "API: extract" begin
    dir = make_test_dir()
    hash = tree_hash(dir)
    tarball = Tar.create(dir)
    rm(dir, recursive=true)
    @test hash == Tar.tree_hash(tarball, skip_empty=true)
    @test hash != Tar.tree_hash(tarball, skip_empty=false)

    @testset "without predicate" begin
        arg_readers(tarball) do tar
            # extract(tarball)
            @arg_test tar begin
                dir = Tar.extract(tar)
                check_tree_hash(hash, dir)
            end
            # extract(tarball, dir) — non-existent
            @arg_test tar begin
                dir = tempname()
                Tar.extract(tar, dir)
                check_tree_hash(hash, dir)
            end
            # extract(tarball, dir) — existent, empty
            @arg_test tar begin
                dir = mktempdir()
                Tar.extract(tar, dir)
                check_tree_hash(hash, dir)
            end
            # extract(tarball, dir) — non-directory (error)
            @arg_test tar begin
                file = tempname()
                touch(file)
                @test_throws ErrorException Tar.extract(tar, file)
                rm(file)
            end
            # extract(tarball, dir) — non-empty directory (error)
            @arg_test tar begin
                dir = mktempdir()
                touch(joinpath(dir, "file"))
                @test_throws ErrorException Tar.extract(tar, dir)
                rm(dir, recursive=true)
            end
        end
    end

    NON_STDLIB_TESTS &&
    @testset "inconvenient stream buffering" begin
        # We will try feeding in an adversarial length that used to cause an assertion error
        open(tarball, read=true) do io
            # This will cause an assertion error because we know the padded space beyond the
            # end of the test file content will be larger than 17 bytes, causing the `for`
            # loop to exit early, failing the assertion.
            @test hash == Tar.tree_hash(ChaosBufferStream(io; chunksizes=[17]); skip_empty=true)
        end

        # This also affected read_data()
        mktempdir() do dir
            open(tarball, read=true) do io
                Tar.extract(ChaosBufferStream(io; chunksizes=[17]), dir)
                check_tree_hash(hash, dir)
            end
        end

        # We also perform a fuzzing test to convince ourselves there are no other errors
        # of this type within `Tar.tree_hash()`.
        for idx in 1:100
            open(tarball, read=true) do io
                @test hash == Tar.tree_hash(ChaosBufferStream(io), skip_empty=true)
            end
        end
    end

    @testset "with predicate" begin
        # generate a version of dir with .skip entries
        dir = make_test_dir(true)
        tarball = Tar.create(dir)
        rm(dir, recursive=true)
        @test hash != Tar.tree_hash(tarball, skip_empty=true)
        @test hash != Tar.tree_hash(tarball, skip_empty=false)

        # predicate to skip paths ending in `.skip`
        predicate = hdr -> !any(splitext(p)[2] == ".skip" for p in split(hdr.path, '/'))
        @test hash == Tar.tree_hash(predicate, tarball, skip_empty=true)
        @test hash != Tar.tree_hash(predicate, tarball, skip_empty=false)

        arg_readers(tarball) do tar
            # extract(predicate, tarball)
            @arg_test tar begin
                dir = Tar.extract(predicate, tar)
                check_tree_hash(hash, dir)
            end
            # extract(predicate, tarball, dir) — non-existent
            @arg_test tar begin
                dir = tempname()
                Tar.extract(predicate, tar, dir)
                check_tree_hash(hash, dir)
            end
            # extract(predicate, tarball, dir) — existent, empty
            @arg_test tar begin
                dir = mktempdir()
                Tar.extract(predicate, tar, dir)
                check_tree_hash(hash, dir)
            end
            # extract(predicate, tarball, dir) — non-directory (error)
            @arg_test tar begin
                file = tempname()
                touch(file)
                @test_throws ErrorException Tar.extract(predicate, tar, file)
                rm(file)
            end
            # extract(predicate, tarball, dir) — non-empty directory (error)
            @arg_test tar begin
                dir = mktempdir()
                touch(joinpath(dir, "file"))
                @test_throws ErrorException Tar.extract(predicate, tar, dir)
                rm(dir, recursive=true)
            end
        end
    end
end

@testset "API: rewrite" begin
    # reference standard tarball
    reference, hash₁ = make_test_tarball()
    ref = read(reference)

    # alternate format tarball
    if @isdefined(gtar)
        # alternate tarball made by GNU tar
        alternate, hash₂ = make_test_tarball() do root
            tarball = tempname()
            gtar(gtar -> run(`$gtar -C $root -cf $tarball .`))
            return tarball
        end
        @test hash₁ == hash₂
        @test ref != read(alternate)
    else
        # at least test the plumbing
        alternate = tempname()
        cp(reference, alternate)
    end

    @testset "without predicate" begin
        for tarball in (reference, alternate)
            arg_readers(tarball) do old
                # rewrite(old)
                @arg_test old begin
                    new_file = Tar.rewrite(old)
                    @test ref == read(new_file)
                    rm(new_file)
                end
                # rewrite(old, new)
                arg_writers() do new_file, new
                    @arg_test old new begin
                        @test new == Tar.rewrite(old, new)
                    end
                    @test ref == read(new_file)
                end
            end
        end
    end

    @testset "with predicate" begin
        # made up order-independent tarball predicate
        predicate = hdr ->
            hdr.type == :symlink ? isodd(length(hdr.link)) : isodd(hdr.size)
        filtered = Tar.create(Tar.extract(predicate, reference))
        ref = read(filtered)
        rm(filtered)

        for tarball in (reference, alternate)
            arg_readers(tarball) do old
                # rewrite(predicate, old)
                @arg_test old begin
                    new_file = Tar.rewrite(predicate, old)
                    @test ref == read(new_file)
                    rm(new_file)
                end
                # rewrite(predicate, old, new)
                arg_writers() do new_file, new
                    @arg_test old new begin
                        @test new == Tar.rewrite(predicate, old, new)
                    end
                    @test ref == read(new_file)
                end
            end
        end
    end

    # cleanup
    rm(alternate)
    rm(reference)
end

@testset "API: skeletons" begin
    # make some tarballs to test with
    tarballs = Dict{String,Bool}() # value indicates if we generated
    let dir = make_test_dir()
        tarballs[Tar.create(dir)] = true
        rm(dir, recursive=true)
    end
    tarballs[make_test_tarball()[1]] = true
    if @isdefined(gtar)
        tarball, _ = make_test_tarball() do root
            tarball = tempname()
            gtar(gtar -> run(`$gtar --format=gnu -C $root -cf $tarball .`))
            return tarball
        end
        tarballs[tarball] = false # not generated by Tar.jl
    end
    for (tarball, flag) in collect(tarballs)
        tarball′ = tempname()
        cp(tarball, tarball′)
        open(tarball′, append=true) do io
            write(io, zeros(UInt8, 1024))
            write(io, rand(UInt8, 666))
        end
        tarballs[tarball′] = flag
    end

    for (tarball, flag) in tarballs
        reference = read(tarball)
        # first, generate a skeleton
        skeleton = tempname()
        dir = Tar.extract(tarball, skeleton=skeleton)
        @test isfile(skeleton)
        # test skeleton listing
        hdrs = Tar.list(tarball)
        arg_readers(skeleton) do skel
            @arg_test skel @test hdrs == Tar.list(skel)
        end
        if flag && @isdefined(gtar)
            # GNU tar can list skeleton files of tarballs we generated
            paths = sort!([hdr.path for hdr in hdrs])
            @test paths == sort!(gtar(gtar -> readlines(`$gtar -tf $skeleton`)))
        end
        hdrs = Tar.list(tarball, raw=true)
        arg_readers(skeleton) do skel
            @arg_test skel @test hdrs == Tar.list(skel, raw=true)
            # test reconstruction from skeleton
            @arg_test skel begin
                tarball′ = Tar.create(dir, skeleton=skel)
                @test reference == read(tarball′)
                rm(tarball′)
            end
        end
        # check that extracting skeleton to IO works
        arg_writers() do skeleton′, skel
            @arg_test skel Tar.extract(tarball, skeleton=skel)
            @test read(skeleton) == read(skeleton′)
        end
        rm(skeleton)
    end

    # cleanup
    foreach(rm, keys(tarballs))
end

if Sys.iswindows() && Sys.which("icacls") !== nothing && VERSION >= v"1.6"
    @testset "windows permissions" begin
        tarball, hash = make_test_tarball()
        mktempdir() do dir
            Tar.extract(tarball, dir)
            f_path = joinpath(dir, "0-ffffffff")
            @test isfile(f_path)
            @test !Sys.isexecutable(f_path)

            x_path = joinpath(dir, "0-xxxxxxxx")
            @test isfile(x_path)
            @test Sys.isexecutable(x_path)

            f_acl = readchomp(`icacls $(f_path)`)
            @test occursin("Everyone:(R,WA)", f_acl)
            x_acl = readchomp(`icacls $(x_path)`)
            @test occursin("Everyone:(RX,WA)", x_acl)
        end
    end
end


text:
padfac(alg::NUFT) = alg.padfac

function padimage(::NUFT, img)
    #pf = padfac(alg)
    #cimg = convert(Matrix{Complex{eltype(img)}}, img.img)
    return IntensityMap(img, img.fovx, img.fovy, img.pulse)
    # ny,nx = size(img)
    # nnx = nextpow(2, pf*nx)
    # nny = nextpow(2, pf*ny)
    # nsx = nnx÷2-nx÷2
    # nsy = nny÷2-ny÷2
    # cimg = convert(Matrix{Complex{eltype(img)}}, img)
    # return PaddedView(zero(eltype(cimg)), cimg,
    #                   (1:nnx, 1:nny),
    #                   (nsx+1:nsx+nx, nsy+1:nsy+ny)
    #                  )
end

padimage(alg::ObservedNUFT, img) = padimage(alg.alg, img)

function create_cache(alg::ObservedNUFT, img)
    pimg = padimage(alg, img)

    # make nuft plan
    dx, dy = pixelsizes(img)
    plan = plan_nuft(alg, pimg, dx, dy)
    # get phases and pulse functions
    phases = make_phases(alg, img)

    return create_cache(alg, plan, phases, pimg)
end

function create_cache(alg::NUFT, img)
    pimg = padimage(alg, img)
    return NUFTCache(alg, nothing, nothing, img.pulse, pimg)
end

function update_cache(cache::NUFTCache, img)
    pimg = padimage(cache.alg, img)
    cache2 = update_phases(cache, img)
    create_cache(cache2.alg, cache2.plan, cache2.phases, pimg)
end

function update_phases(cache::NUFTCache, img)
    if cache.pulse !== img.pulse
        phases = make_phases(cache.alg, img)
        return @set cache.phases = phases
    else
        return cache
    end
end

function nocachevis(m::ModelImage{M,I,<:NUFTCache}, u::AbstractArray, v::AbstractArray) where {M,I}
    alg = ObservedNUFT(m.cache.alg, vcat(u', v'))
    cache = create_cache(alg, m.image)
    m = @set m.cache = cache
    return _visibilities(m, u, v)
end

function checkuv(uv, u, v)
    @assert u == @view(uv[1,:]) "Specified u don't match uv in cache. Did you pass the correct u,v?"
    @assert v == @view(uv[2,:]) "Specified v don't match uv in cache. Did you pass the correct u,v?"
end

function nuft(A, b)
    return A*(b .+ 0im)
end

# function ChainRulesCore.rrule(::typeof(nuft), A, b)
#     vis = A*b
#     println("Pld")
#     function nuft_pullback(Δy)
#         Δf = NoTangent()
#         ΔA = @thunk(A'*Δy)
#         Δb = @thunk(Δy*b')
#         return Δf, ΔA, Δb
#     end
#     return vis, nuft_pullback
# end

# ReverseDiff.@grad_from_chainrules nuft(A::ReverseDiff.TrackedArray, b::ReverseDiff.TrackedArray)
# ReverseDiff.@grad_from_chainrules nuft(A, b::ReverseDiff.TrackedArray)
# ReverseDiff.@grad_from_chainrules nuft(A::ReverseDiff.TrackedArray, b)
# ReverseDiff.@grad_from_chainrules nuft(A, b::Vector{<:ReverseDiff.TrackedReal})

function _visibilities(m::ModelImage{M,I,<:NUFTCache{A}},
                      u::AbstractArray,
                      v::AbstractArray) where {M,I,A<:ObservedNUFT}
    checkuv(m.cache.alg.uv, u, v)
    #vr = real.(m.cache.plan)
    #vi = imag.(m.cache.plan)
    vis =  nuft(m.cache.plan, m.cache.img)
    conj.(vis).*m.cache.phases
    #return vis
end

# function _visibilities(m::ModelImage{M,I,<:NUFTCache{A}},
#                        u::AbstractArray,
#                        v::AbstractArray) where {M,I,G<:DFTAlg,A<:ObservedNUFT{G}}
#     checkuv(m.cache.alg.uv, u, v)
#     #vr = real.(m.cache.plan)
#     #vi = imag.(m.cache.plan)
#     visr = nuft(m.cache.plan.re, m.cache.img)
#     visi = nuft(m.cache.plan.im, m.cache.img)
#     vis =  StructArray{Complex{eltype(visr)}}((visr, visi))
#     #vis = m.cache.plan*m.cache.img
#     conj.(vis).*m.cache.phases
#     #return vis
# end

function _visibilities(m::ModelImage{M,I,<:NUFTCache{A}},
                      u::AbstractArray,
                      v::AbstractArray) where {M,I,A<:NUFT}
    return nocachevis(m, u, v)
end

function _amplitudes(m::ModelImage{M,I,C}, u::AbstractArray, v::AbstractArray) where {M,I,C<:NUFTCache}
    vis = visibilities(m, u, v)
    return abs.(vis)
end

# function _bispectra(m::ModelImage{M,I,C},
#                     u1::AbstractArray,
#                     v1::AbstractArray,
#                     u2::AbstractArray,
#                     v2::AbstractArray,
#                     u3::AbstractArray,
#                     v3::AbstractArray) where {M,I,C<:NUFTCache}
#     # TODO: Fix this so we can actually use a stored cache for closures
#     n = length(u1)
#     # In testing it is faster to concatenate everything into 1 vector and then
#     # take an NFFT of that
#     u = vcat(u1,u2,u3)
#     v = vcat(v1,v2,v3)
#     vis = nocachevis(m, u, v)
#     @inbounds (@view(vis[1:n])).*(@view(vis[n+1:2n])).*(@view vis[2n+1:3n])
# end

# function _closure_phases(m::ModelImage{M,I,C},
#                         u1::AbstractArray,
#                         v1::AbstractArray,
#                         u2::AbstractArray,
#                         v2::AbstractArray,
#                         u3::AbstractArray,
#                         v3::AbstractArray) where {M,I,C<:NUFTCache}
#     # TODO: Fix this so we can actually use a stored cache for closures
#     n = length(u1)
#     # In testing it is faster to concatenate everything into 1 vector and then
#     # take an NFFT of that
#     u = vcat(u1,u2,u3)
#     v = vcat(v1,v2,v3)
#     vis = nocachevis(m, u, v)
#     phase = @inbounds angle.(@view(vis[1:n])) .+ angle.(@view(vis[n+1:2n])) .+ angle.(@view vis[2n+1:3n])
#     return phase
# end

# function _logclosure_amplitudes(m::ModelImage{M,I,C},
#                                u1::AbstractArray,
#                                v1::AbstractArray,
#                                u2::AbstractArray,
#                                v2::AbstractArray,
#                                u3::AbstractArray,
#                                v3::AbstractArray,
#                                u4::AbstractArray,
#                                v4::AbstractArray,
#                                ) where {M,I,C<:NUFTCache}
#     # TODO: Fix this so we can actually use a stored cache for closures
#     n = length(u1)
#     u = vcat(u1,u2,u3,u4)
#     v = vcat(v1,v2,v3,v4)
#     vis = nocachevis(m, u, v)
#     amp1 = @view vis[1:n]
#     amp2 = @view vis[n+1:2n]
#     amp3 = @view vis[2n+1:3n]
#     amp4 = @view vis[3n+1:4n]
#     lcamp = @. log(abs(amp1)*abs(amp2)/(abs(amp3)*abs(amp4)))
#     return lcamp
# end

"""
    NFFTAlg
Uses a non-uniform FFT to compute the visibilities.
You can optionally pass uv which are the uv positions you will
compute the NFFT at. This can allow for the NFFT plan to be cached improving
performance

# Fields
$(FIELDS)

"""
Base.@kwdef struct NFFTAlg <: NUFT
    """
    Amount to pad the image
    """
    padfac::Int = 1
    """
    Controls the accuracy of the NFFT usually don't need to change this
    """
    m::Int = 10
end
include(joinpath(@__DIR__, "nfft_alg.jl"))

"""
    DFTAlg
Uses a discrete fourier transform. This is not very efficient for larger images. In those cases
 NFFTAlg or FFTAlg are more reasonable. For small images this is a reasonable choice especially
since it's easy to define derivatives.
"""
struct DFTAlg <: NUFT end
include(joinpath(@__DIR__, "dft_alg.jl"))


text:
using NiLang, NiLang.AD

@i function isqrt(out!, x::T; atol::Float64=1e-8) where T
    @anc anc1::T
    @anc anc2::T
    @anc yn::T
    @anc ynpre::T
    @anc diff::T
    @anc count::Int

    @routine forward begin
        yn += x
        while (!isapprox(val(yn)^2, val(x); atol=atol) && count<50, count!=0)
            anc1 += x / yn
            anc1 -= yn
            yn += anc1 / 2
            # speudo inverse
            @routine getynpre begin
                @safe @show anc2
                if (val(x) > 1.0, ~)
                    anc2 += yn^2
                    anc2 -= x
                    ynpre += anc2^0.5  # +-
                else
                    anc2 -= yn^2
                    anc2 += x
                    ynpre -= anc2^0.5  # +-
                end
                ynpre += yn
                @safe @show ynpre, yn, anc2
            end
            @safe @show ynpre
            anc1 += ynpre
            anc1 -= x/ynpre
            ## uncompute ynpre and anc2
            ~@routine getynpre
            count += 1
        end
    end
    out! += yn

    ~@routine forward
end

using Test
# NOTE: to allow high performance use of f += T(x).
# Now, this kind of use is a performance killer.
@testset "isqrt" begin
    out = 0.0
    x = 1.3
    @instr isqrt(out, x)
    res = sqrt(x)
    @test check_inv(isqrt,(out, x))
    @test isapprox(out, res; atol=1e-8)

    @show "@"
    out = 0.0
    x = 0.3
    @instr isqrt(out, x)
    res = sqrt(x)
    @test check_inv(isqrt,(out, x); verbose=true)
    @test isapprox(out, res; atol=1e-8)

    out = 0.0
    x = 1.0

    @instr isqrt(out, x)
    res = sqrt(x)
    @test check_inv(isqrt,(out, x))
    @test out ≈ res
end

@testset "isqrt grad" begin
    out = 0.0
    x = 1.6
    gres = sqrt(x)
    @test check_inv(isqrt, (Loss(out), x); verbose=true)
    @test check_grad(isqrt, (Loss(out), x); verbose=true)

    out = 0.0
    x = 1.6
    @instr isqrt'(Loss(out), x)
    @test grad(x) ≈ gres
end