text:
"""
    RBM{V,H,W}

RBM, with visible layer of type `V`, hidden layer of type `H`, and weights of type `W`.
"""
struct RBM{V,H,W}
    visible::V
    hidden::H
    w::W
    """
        RBM(visible, hidden, w)

    Creates a Restricted Boltzmann machine with `visible` and `hidden` layers and weights `w`.
    """
    function RBM(visible::AbstractLayer, hidden::AbstractLayer, w::AbstractArray)
        @assert size(w) == (size(visible)..., size(hidden)...)
        return new{typeof(visible), typeof(hidden), typeof(w)}(visible, hidden, w)
    end
end

visible(rbm::RBM) = rbm.visible
hidden(rbm::RBM) = rbm.hidden
weights(rbm::RBM) = rbm.w
flat_w(rbm::RBM) = reshape(weights(rbm), length(visible(rbm)), length(hidden(rbm)))
flat_v(rbm::RBM, v::AbstractArray) = flatten(visible(rbm), v)
flat_h(rbm::RBM, h::AbstractArray) = flatten(hidden(rbm), h)

"""
    inputs_h_from_v(rbm, v)

Interaction inputs from visible to hidden layer.
"""
function inputs_h_from_v(rbm::RBM, v::AbstractArray)
    wflat = flat_w(rbm)
    vflat = activations_convert_maybe(wflat, flat_v(rbm, v))
    iflat = wflat' * vflat
    return reshape(iflat, size(hidden(rbm))..., batch_size(visible(rbm), v)...)
end

"""
    inputs_v_from_h(rbm, h)

Interaction inputs from hidden to visible layer.
"""
function inputs_v_from_h(rbm::RBM, h::AbstractArray)
    wflat = flat_w(rbm)
    hflat = activations_convert_maybe(wflat, flat_h(rbm, h))
    iflat = wflat * hflat
    return reshape(iflat, size(visible(rbm))..., batch_size(hidden(rbm), h)...)
end

# compat
inputs_v_to_h(rbm::RBM, v::AbstractArray) = inputs_h_from_v(rbm, v)
inputs_h_to_v(rbm::RBM, h::AbstractArray) = inputs_v_from_h(rbm, h)

"""
    free_energy(rbm, v)

Free energy of visible configuration (after marginalizing hidden configurations).
"""
function free_energy(rbm::RBM, v::AbstractArray)
    E = energy(visible(rbm), v)
    inputs = inputs_h_from_v(rbm, v)
    F = free_energy(hidden(rbm), inputs)
    return E + F
end

"""
    energy(rbm, v, h)

Energy of the rbm in the configuration `(v,h)`.
"""
function energy(rbm::RBM, v::AbstractArray, h::AbstractArray)
    Ev = energy(visible(rbm), v)
    Eh = energy(hidden(rbm), h)
    Ew = interaction_energy(rbm, v, h)
    return Ev .+ Eh .+ Ew
end

"""
    interaction_energy(rbm, v, h)

Weight mediated interaction energy.
"""
function interaction_energy(rbm::RBM, v::AbstractArray, h::AbstractArray)
    bsz = batch_size(rbm, v, h)
    if ndims(visible(rbm)) == ndims(v) || ndims(hidden(rbm)) == ndims(h)
        E = -flat_v(rbm, v)' * flat_w(rbm) * flat_h(rbm, h)
    elseif length(visible(rbm)) ≥ length(hidden(rbm))
        inputs = inputs_h_from_v(rbm, v)
        E = -sum(inputs .* h; dims = 1:ndims(hidden(rbm)))
    else
        inputs = inputs_v_from_h(rbm, h)
        E = -sum(v .* inputs; dims=1:ndims(visible(rbm)))
    end
    return reshape_maybe(E, bsz)
end

"""
    sample_h_from_v(rbm, v)

Samples a hidden configuration conditional on the visible configuration `v`.
"""
function sample_h_from_v(rbm::RBM, v::AbstractArray)
    inputs = inputs_h_from_v(rbm, v)
    return sample_from_inputs(hidden(rbm), inputs)
end

"""
    sample_v_from_h(rbm, h)

Samples a visible configuration conditional on the hidden configuration `h`.
"""
function sample_v_from_h(rbm::RBM, h::AbstractArray)
    inputs = inputs_v_from_h(rbm, h)
    return sample_from_inputs(visible(rbm), inputs)
end

"""
    sample_v_from_v(rbm, v; steps = 1)

Samples a visible configuration conditional on another visible configuration `v`.
Ensures type stability by requiring that the returned array is of the same type as `v`.
"""
function sample_v_from_v(rbm::RBM, v::AbstractArray; steps::Int = 1)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    for _ in 1:steps
        v = oftype(v, sample_v_from_v_once(rbm, v))
    end
    return v
end

"""
    sample_h_from_h(rbm, h; steps = 1)

Samples a hidden configuration conditional on another hidden configuration `h`.
Ensures type stability by requiring that the returned array is of the same type as `h`.
"""
function sample_h_from_h(rbm::RBM, h::AbstractArray; steps::Int = 1)
    @assert size(hidden(rbm)) == size(h)[1:ndims(hidden(rbm))]
    for _ in 1:steps
        h = oftype(h, sample_h_from_h_once(rbm, h))
    end
    return h
end

function sample_v_from_v_once(rbm::RBM, v::AbstractArray)
    h = sample_h_from_v(rbm, v)
    v = sample_v_from_h(rbm, h)
    return v
end

function sample_h_from_h_once(rbm::RBM, h::AbstractArray)
    v = sample_v_from_h(rbm, h)
    h = sample_h_from_v(rbm, v)
    return h
end

"""
    mean_h_from_v(rbm, v)

Mean unit activation values, conditioned on the other layer, <h | v>.
"""
function mean_h_from_v(rbm::RBM, v::AbstractArray)
    inputs = inputs_h_from_v(rbm, v)
    return mean_from_inputs(hidden(rbm), inputs)
end

"""
    mean_v_from_h(rbm, v)

Mean unit activation values, conditioned on the other layer, <v | h>.
"""
function mean_v_from_h(rbm::RBM, h::AbstractArray)
    inputs = inputs_v_from_h(rbm, h)
    return mean_from_inputs(visible(rbm), inputs)
end

"""
    var_v_from_h(rbm, v)

Variance of unit activation values, conditioned on the other layer, var(v | h).
"""
function var_v_from_h(rbm::RBM, h::AbstractArray)
    inputs = inputs_v_from_h(rbm, h)
    return var_from_inputs(visible(rbm), inputs)
end

"""
    var_h_from_v(rbm, v)

Variance of unit activation values, conditioned on the other layer, var(h | v).
"""
function var_h_from_v(rbm::RBM, v::AbstractArray)
    inputs = inputs_h_from_v(rbm, v)
    return var_from_inputs(hidden(rbm), inputs)
end

"""
    mode_v_from_h(rbm, h)

Mode unit activations, conditioned on the other layer.
"""
function mode_v_from_h(rbm::RBM, h::AbstractArray)
    inputs = inputs_v_from_h(rbm, h)
    return mode_from_inputs(visible(rbm), inputs)
end

"""
    mode_h_from_v(rbm, v)

Mode unit activations, conditioned on the other layer.
"""
function mode_h_from_v(rbm::RBM, v::AbstractArray)
    inputs = inputs_h_from_v(rbm, v)
    return mode_from_inputs(hidden(rbm), inputs)
end

"""
    batch_size(rbm, v, h)

Returns the batch size if `energy(rbm, v, h)` were computed.
"""
function batch_size(rbm::RBM, v::AbstractArray, h::AbstractArray)
    v_bsz = batch_size(visible(rbm), v)
    h_bsz = batch_size(hidden(rbm), h)
    if isempty(v_bsz)
        return h_bsz
    elseif isempty(h_bsz)
        return v_bsz
    else
        return join_batch_size(v_bsz, h_bsz)
    end
end

function join_batch_size(bsz_1::Tuple{Int,Vararg{Int}}, bsz_2::Tuple{Int,Vararg{Int}})
    if length(bsz_1) > length(bsz_2)
        D = length(bsz_2)
        sz2 = bsz_1[(D + 1):end]
    else
        D = length(bsz_1)
        sz2 = bsz_2[(D + 1):end]
    end
    sz1 = map(bsz_1[1:D], bsz_2[1:D]) do b1, b2
        bmin, bmax = minmax(b1, b2)
        @assert bmin == 1 || bmin == bmax
        bmax
    end
    return (sz1..., sz2...)
end

"""
    reconstruction_error(rbm, v; steps = 1)

Stochastic reconstruction error of `v`.
"""
function reconstruction_error(rbm::RBM, v::AbstractArray; steps::Int=1)
    @assert size(visible(rbm)) == size(v)[1:ndims(visible(rbm))]
    v1 = sample_v_from_v(rbm, v; steps)
    ϵ = mean(abs.(v .- v1); dims = 1:ndims(visible(rbm)))
    if ndims(v) == ndims(visible(rbm))
        return only(ϵ)
    else
        return reshape(ϵ, batch_size(visible(rbm), v))
    end
end

"""
    mirror(rbm)

Returns a new RBM with visible and hidden layers flipped.
"""
function mirror(rbm::RBM)
    p(i::Int) = i ≤ ndims(visible(rbm)) ? i + ndims(hidden(rbm)) : i - ndims(visible(rbm))
    perm = ntuple(p, ndims(weights(rbm)))
    w = permutedims(weights(rbm), perm)
    return RBM(hidden(rbm), visible(rbm), w)
end

"""
    ∂free_energy(rbm, v)

Gradient of `free_energy(rbm, v)` with respect to model parameters.
If `v` consists of multiple samples (batches), then an average is taken.
"""
function ∂free_energy(
    rbm::RBM, v::AbstractArray; wts = nothing, stats = suffstats(rbm, v; wts)
)
    inputs = inputs_h_from_v(rbm, v)
    ∂v = ∂energy(visible(rbm), stats)
    ∂Γ = ∂free_energies(hidden(rbm), inputs)
    ∂h = map(∂Γ) do ∂f
        batchmean(hidden(rbm), ∂f; wts)
    end
    h = grad2ave(hidden(rbm), ∂Γ)
    ∂w = ∂interaction_energy(rbm, v, h; wts)
    return (visible = ∂v, hidden = ∂h, w = ∂w)
end

function ∂interaction_energy(rbm::RBM, v::AbstractArray, h::AbstractArray; wts = nothing)
    bsz = batch_size(rbm, v, h)
    if ndims(visible(rbm)) == ndims(v) && ndims(hidden(rbm)) == ndims(h)
        wts::Nothing
        ∂wflat = -vec(v) * vec(h)'
    elseif ndims(visible(rbm)) == ndims(v)
        ∂wflat = -vec(v) * vec(batchmean(hidden(rbm), h; wts))'
    elseif ndims(hidden(rbm)) == ndims(h)
        ∂wflat = -vec(batchmean(visible(rbm), v; wts)) * vec(h)'
    else
        hflat = flatten(hidden(rbm), h)
        vflat = activations_convert_maybe(hflat, flatten(visible(rbm), v))
        @assert isnothing(wts) || size(wts) == batch_size(visible(rbm), v)
        if isnothing(wts)
            ∂wflat = -vflat * hflat' / size(vflat, 2)
        else
            @assert size(wts) == bsz
            @assert batch_size(visible(rbm), v) == batch_size(hidden(rbm), h) == size(wts)
            ∂wflat = -vflat * Diagonal(vec(wts)) * hflat' / sum(wts)
        end
    end
    ∂w = reshape(∂wflat, size(weights(rbm)))
    return ∂w
end

suffstats(rbm::RBM, v::AbstractArray; wts=nothing) = suffstats(visible(rbm), v; wts)


text:
using Documenter, Tensors

# Setup for doctests in docstrings
DocMeta.setdocmeta!(Tensors, :DocTestSetup,
    quote
        using Random
        Random.seed!(1234)
        using Tensors
    end
)

makedocs(
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    modules = [Tensors],
    sitename = "Tensors.jl",
    pages = Any[
        "Home" => "index.md",
        "Manual" => [
            "man/constructing_tensors.md",
            "man/indexing.md",
            "man/binary_operators.md",
            "man/other_operators.md",
            "man/storing_tensors.md",
            "man/automatic_differentiation.md",
            ],
        "Benchmarks" => "benchmarks.md",
        "Demos" => "demos.md"
    ]
)

deploydocs(
    repo = "github.com/KristofferC/Tensors.jl.git",
)

ext: jl


text:
module MoveBlock

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict()

buttonColor = (71, 64, 112, 255) ./ 255

directions = Maple.move_block_directions
for direction in directions, steerable in false:true, fast in false:true
    key = "Move Block ($(uppercasefirst(direction)), $(fast ? "Fast" : "Slow")$(steerable ? ", Steerable" :    ""))"
    placements[key] = Ahorn.EntityPlacement(
        Maple.MoveBlock,
        "rectangle",
        Dict{String, Any}(
            "canSteer" => steerable,
            "direction" => direction,
            "fast" => fast
        )
    )
end

Ahorn.editingOptions(entity::Maple.MoveBlock) = Dict{String, Any}(
    "direction" => Maple.move_block_directions
)
Ahorn.minimumSize(entity::Maple.MoveBlock) = 16, 16
Ahorn.resizable(entity::Maple.MoveBlock) = true, true

Ahorn.selection(entity::Maple.MoveBlock) = Ahorn.getEntityRectangle(entity)

midColor = (4, 3, 23) ./ 255
highlightColor = (59, 50, 101) ./ 255

arrows = Dict{String, String}(
    "up" => "objects/moveBlock/arrow02",
    "left" => "objects/moveBlock/arrow04",
    "right" => "objects/moveBlock/arrow00",
    "down" => "objects/moveBlock/arrow06",
)

button = "objects/moveBlock/button"

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.MoveBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    canSteer = get(entity.data, "canSteer", false)
    direction = lowercase(get(entity.data, "direction", "up"))
    arrowSprite = Ahorn.getSprite(arrows[lowercase(direction)], "Gameplay")

    frame = "objects/moveBlock/base"
    if canSteer
        if direction == "up" || direction == "down"
            frame = "objects/moveBlock/base_v"

        else
            frame = "objects/moveBlock/base_h"
        end
    end

    Ahorn.drawRectangle(ctx, 2, 2, width - 4, height - 4, highlightColor, highlightColor)
    Ahorn.drawRectangle(ctx, 8, 8, width - 16, height - 16, midColor)

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, 0, 8, 0, 8, 8)
        Ahorn.drawImage(ctx, frame, (i - 1) * 8, height - 8, 8, 16, 8, 8)

        if canSteer && (direction != "up" && direction != "down")
            Ahorn.drawImage(ctx, button, (i - 1) * 8, -2, 6, 0, 8, 6, tint=buttonColor)
        end
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, frame, 0, (i - 1) * 8, 0, 8, 8, 8)
        Ahorn.drawImage(ctx, frame, width - 8, (i - 1) * 8, 16, 8, 8, 8)

        if canSteer && (direction == "up" || direction == "down")
            Ahorn.Cairo.save(ctx)

            Ahorn.rotate(ctx, -pi / 2)
            Ahorn.drawImage(ctx, button, i * 8 - height - 8, -2, 6, 0, 8, 6, tint=buttonColor)
            Ahorn.scale(ctx, 1, -1)
            Ahorn.drawImage(ctx, button, i * 8 - height - 8, -2 - width, 6, 0, 8, 6, tint=buttonColor)

            Ahorn.Cairo.restore(ctx)
        end
    end

    Ahorn.drawImage(ctx, frame, 0, 0, 0, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, 0, 16, 0, 8, 8)
    Ahorn.drawImage(ctx, frame, 0, height - 8, 0, 16, 8, 8)
    Ahorn.drawImage(ctx, frame, width - 8, height - 8, 16, 16, 8, 8)

    if canSteer && (direction != "up" && direction != "down")
        Ahorn.Cairo.save(ctx)

        Ahorn.drawImage(ctx, button, 2, -2, 0, 0, 6, 6, tint=buttonColor)
        Ahorn.scale(ctx, -1, 1)
        Ahorn.drawImage(ctx, button, 2 - width, -2, 0, 0, 6, 6, tint=buttonColor)

        Ahorn.Cairo.restore(ctx)
    end

    if canSteer && (direction == "up" || direction == "down")
        Ahorn.Cairo.save(ctx)

        Ahorn.rotate(ctx, -pi / 2)
        Ahorn.drawImage(ctx, button, -height + 2, -2, 0, 0, 8, 6, tint=buttonColor)
        Ahorn.drawImage(ctx, button, -10, -2, 14, 0, 8, 6, tint=buttonColor)
        Ahorn.scale(ctx, 1, -1)
        Ahorn.drawImage(ctx, button, -height + 2, -2 -width, 0, 0, 8, 6, tint=buttonColor)
        Ahorn.drawImage(ctx, button, -10, -2 -width, 14, 0, 8, 6, tint=buttonColor)

        Ahorn.Cairo.restore(ctx)
    end

    Ahorn.drawRectangle(ctx, div(width - arrowSprite.width, 2) + 1, div(height - arrowSprite.height, 2) + 1, 8, 8, highlightColor, highlightColor)
    Ahorn.drawImage(ctx, arrowSprite, div(width - arrowSprite.width, 2), div(height - arrowSprite.height, 2))
end

end


text:
"""
Representation of a geospatial coordinates.

- `lat::Float64`: Latitude.
- `lon::Float64`: Longitude.
- `alt::Float64`: Altitude.
"""
@with_kw struct GeoLocation
    lat::Float64
    lon::Float64
    alt::Float64 = 0.0
end

GeoLocation(lat::AbstractFloat, lon::AbstractFloat)::GeoLocation = GeoLocation(lat=lat, lon=lon)
GeoLocation(point::Vector{<:AbstractFloat})::GeoLocation = GeoLocation(point...)
GeoLocation(point_vector::Vector{<:Vector{<:AbstractFloat}})::Vector{GeoLocation} = [GeoLocation(p...) for p in point_vector]

function Base.:(==)(loc1::GeoLocation, loc2::GeoLocation)
    return loc1.lat == loc2.lat && loc1.lon == loc2.lon && loc1.alt == loc2.alt
end
function Base.hash(loc::GeoLocation, h::UInt)
    for field in fieldnames(GeoLocation)
		h = hash(getproperty(loc, field), h)
	end
    return h
end

"""
OpenStreetMap node.

`T<:Integer`
- `id::T`: OpenStreetMap node id.
- `nodes::Vector{T}`: Node's GeoLocation.
- `tags::AbstractDict{String,Any}`: Metadata tags.
"""
struct Node{T <: Integer}
    id::T
    location::GeoLocation
    tags::Union{Dict{String,Any},Nothing}
end

"""
OpenStreetMap way.

`T<:Integer`
- `id::T`: OpenStreetMap way id.
- `nodes::Vector{T}`: Ordered list of node ids making up the way.
- `tags::AbstractDict{String,Any}`: Metadata tags.
"""
struct Way{T <: Integer}
    id::T
    nodes::Vector{T}
    tags::Dict{String,Any}
end
Way(id::T, nodes, tags::Dict{String, Any}) where T <: Integer = Way(id, convert(Vector{T}, nodes), tags)

"""
OpenStreetMap turn restriction (relation).

`T<:Integer`
- `id::T`: OpenStreetMap relation id.
- `type::String`: Either a `via_way` or `via_node` turn restriction.
- `tags::AbstractDict{String,Any}`: Metadata tags.
- `from_way::T`: Incoming way id to the turn restriction.
- `to_way::T`: Outgoing way id to the turn restriction.
- `via_node::Union{T,Nothing}`: Node id at the centre of the turn restriction.
- `via_way::Union{Vector{T},Nothing}`: Way id at the centre of the turn restriction.
- `is_exclusion::Bool`: Turn restrictions such as `no_left_turn`, `no_right_turn` or `no_u_turn`.
- `is_exclusive::Bool`: Turn restrictions such as `striaght_on_only`, `left_turn_only`, `right_turn_only`.
"""
@with_kw struct Restriction{T <: Integer}
    id::T
    type::String
    tags::Dict{String,Any}
    from_way::T
    to_way::T
    via_node::Union{T,Nothing} = nothing
    via_way::Union{Vector{T},Nothing} = nothing
    is_exclusion::Bool = false
    is_exclusive::Bool = false
end

"""
Container for storing OpenStreetMap node, way, relation and graph related obejcts.

`U <: Integer,T <: Integer,W <: Real`
- `nodes::Dict{T,Node{T}}`: Mapping of node ids to node objects.
- `node_coordinates::Vector{Vector{W}}`: Vector of node coordinates [[lat, lon]...], indexed by graph vertices.
- `ways::Dict{T,Way{T}}`: Mapping of way ids to way objects. Previously called `highways`.
- `node_to_index::OrderedDict{T,U}`: Mapping of node ids to graph vertices.
- `index_to_node::OrderedDict{U,T}`: Mapping of graph vertices to node ids.
- `node_to_way::Dict{T,Vector{T}}`: Mapping of node ids to vector of way ids. Previously called `node_to_highway`.
- `edge_to_way::Dict{Vector{T},T}`: Mapping of edges (adjacent node pairs) to way ids. Previously called `edge_to_highway`.
- `restrictions::Dict{T,Restriction{T}}`: Mapping of relation ids to restriction objects.
- `indexed_restrictions::Union{DefaultDict{U,Vector{MutableLinkedList{U}}},Nothing}`: Mapping of via node ids to ordered sequences of restricted node ids.
- `graph::Union{AbstractGraph,Nothing}`: Either DiGraph, StaticDiGraph, SimpleWeightedDiGraph or MetaDiGraph.
- `weights::Union{SparseMatrixCSC{W,U},Nothing}`: Sparse adjacency matrix (weights between graph vertices), either `:distance` (km), `:time` (hours) or `:lane_efficiency` (time scaled by number of lanes).
- `dijkstra_states::Vector{Vector{U}}`: Vector of dijkstra parent states indexed by source vertices, used to retrieve shortest path from source vertex to any other vertex.
- `kdtree::Union{KDTree,Nothing}`: KDTree used to calculate nearest nodes.
- `weight_type::Union{Symbol,Nothing}`: Either `:distance`, `:time` or `:lane_efficiency`.
"""
@with_kw mutable struct OSMGraph{U <: Integer,T <: Integer,W <: Real}
    nodes::Dict{T,Node{T}} = Dict{T,Node{T}}()
    node_coordinates::Vector{Vector{W}} = Vector{Vector{W}}() # needed for astar heuristic
    ways::Dict{T,Way{T}} = Dict{T,Way{T}}()
    node_to_index::OrderedDict{T,U} = OrderedDict{T,U}()
    index_to_node::OrderedDict{U,T} = OrderedDict{U,T}()
    node_to_way::Dict{T,Vector{T}} = Dict{T,Vector{T}}()
    edge_to_way::Dict{Vector{T},T} = Dict{Vector{T},T}()
    restrictions::Dict{T,Restriction{T}} = Dict{T,Restriction{T}}()
    indexed_restrictions::Union{DefaultDict{U,Vector{MutableLinkedList{U}}},Nothing} = nothing
    graph::Union{AbstractGraph,Nothing} = nothing
    weights::Union{SparseMatrixCSC{W,U},Nothing} = nothing
    dijkstra_states::Union{Vector{Vector{U}},Nothing} = nothing
    kdtree::Union{KDTree,Nothing} = nothing
    weight_type::Union{Symbol,Nothing} = nothing
end

function Base.getproperty(g::OSMGraph, field::Symbol)
    # Ensure renaming of "highways" to "ways" is backwards compatible
    if field === :highways
        Base.depwarn("`highways` field is deprecated, use `ways` field instead", :getproperty)
        return getfield(g, :ways)
    elseif field === :node_to_highway
        Base.depwarn("`node_to_highway` field is deprecated, use `node_to_way` field instead", :getproperty)
        return getfield(g, :node_to_way)
    elseif field === :edge_to_highway
        Base.depwarn("`edge_to_highway` field is deprecated, use `edge_to_way` field instead", :getproperty)
        return getfield(g, :edge_to_way)
    else
        return getfield(g, field)
    end
end

"""
OpenStreetMap building polygon.

`T<:Integer`
- `id::T`: OpenStreetMap building way id.
- `nodes::Vector{T}`: Ordered list of node ids making up the building polyogn.
- `is_outer::Bool`: True if polygon is the outer ring of a multi-polygon.
"""
struct Polygon{T <: Integer}
    id::T
    nodes::Vector{Node{T}}
    is_outer::Bool # or inner
end

"""
OpenStreetMap building.

`T<:Integer`
- `id::T`: OpenStreetMap building way id a simple polygon, relation id if a multi-polygon
- `is_relation::Bool`: True if building is a a multi-polygon / relation.
- `polygons::Vector{Polygon{T}}`: List of building polygons, first is always the outer ring.
- `tags::AbstractDict{String,Any}`: Metadata tags.
"""
struct Building{T <: Integer}
    id::T
    is_relation::Bool # or way
    polygons::Vector{Polygon{T}}
    tags::AbstractDict{String,Any}
end


text:
using Plots

include("functions.jl")

assets_path = "$(@__DIR__)/../assets"

plot_fibo(n) = begin
    (0:(n - 1), fibonacci(n)) |> plot
    title!("Fibonacci Sequence")
    xlabel!("N")
    ylabel!("Fibonacci(N)")
    savefig("$assets_path/fibonacci($n).png")
end

factor_data(n, base, filt=nothing) = begin
    0:(n - 1) |>
    (x -> (collect(x), [factor(i, base) for i ∈ fibonacci(n)])) |>
    zip |>
    collect |>
    (x -> filt === nothing ? x : filter(a -> filt(a[2]), x)) |>
    unzip
end

plot_factor(n, base) = begin
    factor_data(n, base) |> plot
    title!("Factor $base in Fibonacci Numbers")
    xlabel!("N")
    ylabel!("Factor $base in Fibonacci(N)")
    savefig("$assets_path/factor($n,$base).png")
end

plot_factor_except(n, base, except) = begin
    factor_data(n, base, x -> x % except ≠ 0) |> plot
    title!("Factor $base in Fibonacci Numbers except exponent $except")
    xlabel!("N")
    ylabel!("Factor $base in Fibonacci(N) except exponent $except")
    savefig("$assets_path/factor($n,$base)except($except).png")
end

plot_factor_only(n, base, only) = begin
    factor_data(n, base, x -> x % only == 0) |> plot
    title!("Factor $base in Fibonacci Numbers only for exponent $only")
    xlabel!("N")
    ylabel!("Factor $base in Fibonacci(N) only for exponent $only")
    savefig("$assets_path/factor($n,$base)only($only).png")
end

plot_factor(n, base::Vector) = begin
    p = factor_data(n, base[begin]) |> (a, b) -> plot(a, b, label="Base $(base[begin])")
    for i ∈ base[begin + 1:end - 1]
        data = factor_data(n, i)
        plot!(p, data..., label="Base $i")
    end
    factor_data(n, base[end]) |> (a, b) -> plot!(a, b, label="Base $(base[end])")
    title!("Factor $base in Fibonacci Numbers")
    xlabel!("N")
    ylabel!("Multiple factors in Fibonacci(N)")
    bases = map(x -> "$x", base) |> x -> join(x, ",")
    savefig("$assets_path/factor($n,[$bases]).png")
end

function plot_distance(n, base1, base2, e, op)
    (x1, _) = factor_data(n, base1) |>
        zip |>
        collect |>
        k -> filter(d -> d[2] == e, k) |>
        unzip
    (x2, _) = factor_data(n, base2) |>
        zip |>
        collect |>
        k -> filter(d -> d[2] == e, k) |>
        unzip

    (length(x1), length(x2)) |>
    min |>
    (l -> (x1[begin:l], x2[begin:l])) |>
    op |>
    enumerate |>
    unzip |>
    plot

    optitle = if (op == -) "Difference" else "Division" end
    title!("$optitle of factors $base1 and $base2 in Fibonacci(N) for exponent $e")
    xlabel!("N")
    ylabel!("$optitle(N)")
    opname = if (op == -) "minus" else "div" end
    savefig("$assets_path/distance($n,$base1,$base2,$e,$opname).png")
end

max_stars_repo_name: manuelelucchi/Fibonacci

max_issues_repo_name: manuelelucchi/Fibonacci

max_forks_repo_name: manuelelucchi/Fibonacci


text:
print("gail-driver installer")
if VERSION != v"0.6.4"
  print("version is wrong. it must be 0.6.4")
  return
end

gail_root = ENV["GAIL_DRIVER_GITHUB"]

Pkg.add("HDF5")
Pkg.build("HDF5")
using HDF5
Pkg.clone(gail_root * "/Vec.jl")
Pkg.build("Vec.jl")
Pkg.clone("https://github.com/JuliaIO/VideoIO.jl")
Pkg.clone("https://github.com/JuliaGraphics/Cairo.jl")
Pkg.checkout("VideoIO.jl", "mindeps")
Pkg.checkout("Cairo.jl", "app")
Pkg.build("VideoIO.jl")
Pkg.build("Cairo.jl")
Pkg.clone(gail_root * "/AutomotiveDrivingModels.jl")
Pkg.build("AutomotiveDrivingModels.jl")
Pkg.clone(gail_root * "/AutoViz.jl")
Pkg.build("AutoViz.jl")
Pkg.clone(gail_root * "/ForwardNets.jl")
Pkg.build("ForwardNets.jl")
Pkg.clone(gail_root * "/NGSIM.jl")
Pkg.build("NGSIM.jl")
Pkg.clone(gail_root * "/Records.jl")
Pkg.build("Records.jl")
Pkg.clone("https://github.com/tawheeler/AutoDrivers.jl")
using Vec
using AutomotiveDrivingModels
using NGSIM
using AutoViz
using Records
using ForwardNets
using AutoDrivers

max_stars_repo_name: yoshihikokuwahara/gail-driver

max_issues_repo_name: yoshihikokuwahara/gail-driver

max_forks_repo_name: yoshihikokuwahara/gail-driver


text:
"""
    ClaytonCopula{d,T}

Fields:
  - θ::Real - parameter

Constructor

    ClaytonCopula(d, θ)

The [Clayton](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Most_important_Archimedean_copulas) copula in dimension ``d`` is parameterized by ``\\theta \\in [-1,\\infty)`` when ``d=2`` and ``\\theta \\in [0,\\infty)`` if ``d>2``. It is an Archimedean copula with generator :

```math
\\phi(t) = \\left(1+\\mathrm{sign}(\\theta)*t\\right)^{-1\\frac{1}{\\theta}}
```
"""
struct ClaytonCopula{d,T} <: ArchimedeanCopula{d}
    θ::T
end
ClaytonCopula(d,θ) = ((θ >= -1)&(d==2) | (θ>0)) ? ClaytonCopula{d,typeof(θ)}(θ) : @error "Incorect θ value : must be positive if d > 2 and > -1 if d ==2" # this constructor must be implementeD.
ϕ(  C::ClaytonCopula,      t) = (1+sign(C.θ)*t)^(-1/C.θ)
ϕ⁻¹(C::ClaytonCopula,      t) = sign(C.θ)*(t^(-C.θ)-1)
radial_dist(C::ClaytonCopula) = Distributions.Gamma(1/C.θ,1)
τ(C::ClaytonCopula) = C.θ/(C.θ+2)
τ⁻¹(::Type{ClaytonCopula},τ) = 2τ/(1-τ)


text:
module CineFiles

using ImageCore: N0f8, N4f12, Gray
using LRUCache

import Base: eltype, length, size, getindex, firstindex, lastindex, iterate
export CineFile

# using FileIO
# add_format(format"CINE", [0x43, 0x49], ".cine")

struct Time64
    fractions::UInt32
    seconds::UInt32
end

struct CineFileHeader
    Type::UInt16
    HeaderSize::UInt16
    Compression::UInt16
    Version::UInt16
    FirstMovieImage::Int32
    TotalImageCount::UInt32
    FirstImageNum::Int32
    ImageCount::UInt32
    OffImageHeader::UInt32
    OffSetup::UInt32
    OffImageOffsets::UInt32
    TriggerTime::Time64
end

struct BitmapInfoHeader
    Size::UInt32
    Width::Int32
    Height::Int32
    Planes::UInt16
    BitCount::UInt16
    Compression::UInt32
    SizeImage::UInt32
    XPelsPerMeter::Int32
    YPelsPerMeter::Int32
    ClrUsed::UInt32
    ClrImportant::UInt32
    FPS::UInt16
end

readstruct(f, S::T) where {T} = only(read!(f, Vector{S}(undef, 1)))

struct CineHeader{T}
    cine::CineFileHeader
    bitmap::BitmapInfoHeader
    imglocs::Vector{Int64}
    imgoffset::Int
    dt::Vector{Float64}
    tmp::Array{T,2}
end

function CineHeader(fname)
    open(fname) do f
        read(f, UInt16) == UInt(18755) || error(basename(fname), " is not a .cine file")
        seek(f, 0)
        cine = readstruct(f, CineFileHeader)
        seek(f, cine.OffImageHeader)
        bitmap = readstruct(f, BitmapInfoHeader)

        bittype = bitmap.BitCount == 8 ? Gray{N0f8} : Gray{N4f12}
        cine.ImageCount > 0 || error("no images exist in file")

        seek(f, cine.OffImageOffsets)
        imglocs = read!(f, Array{Int64}(undef, cine.ImageCount))

        seekstart(f)
        imgoffset = 1
        while read(f, UInt16) != 1002
            imgoffset += 1
        end

        dt = zeros(cine.ImageCount)
        skip(f, 2)
        for i in eachindex(dt)
            fracstart = read(f, UInt32)
            secstart = read(f, UInt32)
            dt[i] = (secstart - cine.TriggerTime.seconds) + ((fracstart - cine.TriggerTime.fractions)/2^32)
        end

        tmp = Array{bittype, 2}(undef, bitmap.Width, bitmap.Height)
        return CineHeader{bittype}(cine, bitmap, imglocs, imgoffset, dt, tmp)
    end
end

Base.eltype(h::CineHeader{T}) where {T} = T

function readframe!(f::IO, frame, h, frameidx)
    frameidx <= h.cine.ImageCount || error("tried to access nonexistent frame $frameidx")
    seek(f, h.imglocs[frameidx])
    skip(f, read(f, UInt32) - 4)
    read!(f, h.tmp)
    frame .= rotl90(h.tmp)
end

readframe!(filename, frame, h, frameidx) = open(f -> readframe!(f, frame, h, frameidx), filename)
readframe(f, h, frameidx) = readframe!(f, rotl90(h.tmp), h, frameidx)

struct CineFile{T}
    path::String
    header::CineHeader{T}
    data::LRU{Int, Array{T,2}}
end

"""
    CineFile(filepath, cachelimit=0.25)

Load header information and create a frame cache for a Phantom .cine file.
    `cachelimit` sets the maximum size of the frame cache as a fraction of
    your system's free RAM. Indexing and iteration is supported for CineFiles.

"""
function CineFile(filepath, cachelimit=0.25)
    header = CineHeader(filepath)
    framesize = Base.summarysize(header.tmp)
    maxcachedframes = ceil(Int, cachelimit*Sys.free_memory() / framesize)
    data = LRU{Int, Array{eltype(header),2}}(maxsize = maxcachedframes)
    return CineFile(filepath, header, data)
end

Base.length(cf::CineFile) = length(cf.header.dt)
Base.eltype(cf::CineFile) = typeof(cf.header.tmp)
Base.size(cf::CineFile) = (length(cf), size(cf.header.tmp, 2), size(cf.header.tmp, 1))

function Base.getindex(cf::CineFile, idx::Int)
    get!(cf.data, idx) do
        readframe(cf.path, cf.header, idx)
    end
end
Base.firstindex(cf::CineFile) = 1
Base.lastindex(cf::CineFile) = length(cf)
Base.getindex(cf::CineFile, I) = [cf[i] for i in I]

Base.iterate(cf::CineFile, state=1) = state > length(cf) ? nothing : (cf[state], state + 1)

end


text:
module SimplePNGs
using CodecZlib
using UnPack
using FixedPointNumbers
using ColorTypes

function load(name)
    png = SimplePNG(chunkify(read(name)))
    data = extract(png)
    return build!(png, data)
end

function loadpng(name)
    png = SimplePNG(chunkify(read(name)))
end

struct Chunk
    length::UInt32
    name::Vector{UInt8}
    data::Vector{UInt8}
    crc::Vector{UInt8}
end

function chunkify(img::Vector{UInt8})
    res = Chunk[]
    @assert length(img) >= 8
    @assert all(Int.(img[1:8]) .== [137, 80, 78, 71, 13, 10, 26, 10]) # PNG signature
    idx = 9
    while idx < length(img)
        ll = casttoint(img, idx)
        idx += 4
        name = img[idx:idx+3]
        idx += 4
        data = img[idx:idx+ll-1]
        idx += ll
        crc = img[idx:idx+3]
        idx += 4
        chunk = Chunk(ll, name, data, crc)
        push!(res, chunk)
    end

    return res
end

# Cast to int
function casttoint(data, idx)
    @inbounds x = 0 | data[idx]
    x = x << 8
    @inbounds x = x | data[idx + 1]
    x = x << 8
    @inbounds x = x | data[idx + 2]
    x = x << 8
    @inbounds x = x | data[idx + 3]
    return x
end

function casttouint16(data, idx)
    @inbounds x = UInt16(data[idx])
    x = x << 8
    @inbounds x = x | data[idx + 1]

    return x
end

struct SimplePNG
    bitdepth::Int
    colourtype::Int
    compression::Int
    filter::Int
    interlace::Int
    width::Int
    height::Int
    byteshift::Int
    chunks::Vector{Chunk}
end

function SimplePNG(chunks::Vector{Chunk})
    hdr = chunks[1] # Faster search
    data = hdr.data
    width = casttoint(data, 1)
    height = casttoint(data, 5)
    bitdepth = data[9]
    colourtype = data[10]
    compression = data[11]
    filter = data[12]
    interlace = data[13]
    byteshift = bytes(colourtype, bitdepth)

    return SimplePNG(bitdepth, colourtype, compression, filter, interlace, width, height, byteshift, chunks)
end

function name(chunk::Chunk)
    String(Char.(chunk.name))
end

function palette(chunk)
    res = Vector{RGB{N0f8}}(undef, length(chunk) ÷ 3)
    for i in 1:length(res)
        res[i] = pixel(RGB{N0f8}, chunk[3i - 2], chunk[3i - 1], chunk[3i])
    end
    return res
end

function firstdata(png)
    idx = 0
    while idx < length(png.chunks)
        idx += 1
        if name(png.chunks[idx]) == "IDAT"
            return idx
        end
    end

    return idx
end

function extract(png)
    idx = firstdata(png)
    @assert idx > 0
    res = png.chunks[idx].data
    for i in idx+1:length(png.chunks)
        if name(png.chunks[i]) == "IDAT"
            append!(res, png.chunks[i].data)
        end
    end

    data = transcode(ZlibDecompressor, res)
    return data
end

function pixel(::Type{Gray{T}}, c) where T
    return Gray{T}(reinterpret(T, c))
end

function pixel(::Type{GrayA{T}}, c, alpha) where T
    return GrayA{T}(reinterpret(T, c),
                    reinterpret(T, alpha))
end

function pixel(::Type{RGB{T}}, r, g, b) where T
    RGB{T}(reinterpret(T, r),
        reinterpret(T, g),
        reinterpret(T, b))
end

function pixel(::Type{RGBA{T}}, r, g, b, alpha) where T
    RGBA{T}(reinterpret(T, r),
        reinterpret(T, g),
        reinterpret(T, b),
        reinterpret(T, alpha))
end

function a(data, idx, i, j, width, byteshift)
    j == 1 ? 0x00 : data[idx - byteshift]
end

function b(data, idx, i, j, width, byteshift)
    i == 1 ? 0x00 : data[idx - byteshift * width - 1]
end

function c(data, idx, i, j, width, byteshift)
    (i == 1) | (j == 1) ? 0x00 : data[idx - byteshift * width - 1 - byteshift]
end

function avg(data, idx, i, j, width, byteshift)
    a1 = a(data, idx, i, j, width, byteshift)
    b1 = b(data, idx, i, j, width, byteshift)
    return UInt8((UInt16(a1) + UInt16(b1)) >> 1)
end

function paeth(data, idx, i, j, width, byteshift)
    a1 = a(data, idx, i, j, width, byteshift)
    b1 = b(data, idx, i, j, width, byteshift)
    c1 = c(data, idx, i, j, width, byteshift)
    p = Int(a1) + Int(b1) - Int(c1)
    pa = abs(p - a1)
    pb = abs(p - b1)
    pc = abs(p - c1)
    if (pa <= pb) & (pa <= pc)
        return a1
    elseif pb <= pc
        return b1
    else
        return c1
    end
end

function bytes(colourtype, bitdepth)
    if colourtype == 0
        bitdepth == 16 ? 2 : 1
    elseif colourtype == 2
        bitdepth == 16 ? 6 : 3
    elseif colourtype == 3
        1
    elseif colourtype == 4
        bitdepth == 16 ? 4 : 2
    else
        bitdepth == 16 ? 8 : 4
    end
end

function scale_width(width, bpb, bitdepth)
    return if bitdepth < 8
        x, y = divrem(width, bpb)
        x + (y > 0)
    else
        width
    end
end

function level2wh(width, height, level, bpb, bitdepth, width1, height1, scale = true)
    if level == 1
        x, y = divrem(width, 8)
        width1 = scale ? scale_width(x + (y > 0), bpb, bitdepth) : x + (y > 0)
        x, y = divrem(height, 8)
        height1 = x + (y > 0)
    elseif level == 2
        x, y = divrem(width, 8)
        width1 = scale ? scale_width(x + (y > 4), bpb, bitdepth) : x + (y > 4)
    elseif level == 3
        x, y = divrem(width, 4)
        width1 = scale ? scale_width(x + (y > 0), bpb, bitdepth) : x + (y > 0)
        x, y = divrem(height, 8)
        height1 = x + (y > 4)
    elseif level == 4
        x, y = divrem(width, 4)
        width1 = scale ? scale_width(x + (y > 2), bpb, bitdepth) : x + (y > 2)
        x, y = divrem(height, 4)
        height1 = x + (y > 0)
    elseif level == 5
        x, y = divrem(width, 2)
        width1 = scale ? scale_width(x + y, bpb, bitdepth) : x + y
        x, y = divrem(height, 4)
        height1 = x + (y > 2)
    elseif level == 6
        width1 = scale ? scale_width(width ÷ 2, bpb, bitdepth) : width ÷ 2
        x, y = divrem(height, 2)
        height1 = x + y
    elseif level == 7
        width1 = scale ? scale_width(width, bpb, bitdepth) : width
        height1 = height ÷ 2
    end

    return width1, height1
end

function process!(data, idx, width, height, byteshift)
    @inbounds for i in 1:height
        ft = data[idx]
        idx += 1
        for j in 1:width
            for k in 1:byteshift
                if ft == 0x01
                    data[idx] += a(data, idx, i, j, width, byteshift)
                elseif ft == 0x02
                    data[idx] += b(data, idx, i, j, width, byteshift)
                elseif ft == 0x03
                    data[idx] += avg(data, idx, i, j, width, byteshift)
                elseif ft == 0x04
                    data[idx] += paeth(data, idx, i, j, width, byteshift)
                end
                idx += 1
            end
        end
    end

    return idx
end

function process!(png, data)
    @unpack colourtype, bitdepth, width, height, byteshift, interlace = png

    idx = 1
    bpb = 8 ÷ bitdepth
    if interlace == 0
        width = scale_width(width, bpb, bitdepth)
        process!(data, idx, width, height, byteshift)
    else
        height1, width1 = 0, 0
        for level in 1:7
            width1, height1 = level2wh(width, height, level, bpb, bitdepth, width1, height1)
            if (width1 > 0) & (height1 > 0)
                idx = process!(data, idx, width1, height1, byteshift)
            end
        end
    end

    return data
end

function prepare_array(colourtype, bitdepth, width, height)
    if colourtype == 0
        if bitdepth == 16
            T = Gray{N0f16}
        else
            T = Gray{N0f8}
        end
    elseif colourtype == 2
        if bitdepth == 16
            T = RGB{N0f16}
        else
            T = RGB{N0f8}
        end
    elseif colourtype == 3
        T = RGB{N0f8}
    elseif colourtype == 4
        if bitdepth == 8
            T = GrayA{N0f8}
        else
            T = GrayA{N0f16}
        end
    else
        if bitdepth == 8
            T = RGBA{N0f8}
        else
            T = RGBA{N0f16}
        end
    end

    return Array{T}(undef, width, height)
end

function build!(png, data)
    @unpack colourtype, bitdepth, width, height, byteshift, interlace = png
    process!(png, data)
    if colourtype == 3
        idx = findfirst(x -> name(x) == "PLTE", png.chunks)
        chunk = png.chunks[idx]
        plte = palette(chunk.data)
    else
        plte = nothing
    end
    res = prepare_array(colourtype, bitdepth, width, height)
    idx = 1

    if interlace == 0
        build!(res, data, idx, width, height, colourtype, bitdepth, plte, 0)
    else
        bpb = 8 ÷ bitdepth
        width1, height1 = 0, 0
        for level in 1:7
            width1, height1 = level2wh(width, height, level, bpb, bitdepth, width1, height1, false)
            if (width1 > 0) & (height1 > 0)
                idx = build!(res, data, idx, width1, height1, colourtype, bitdepth, plte, level)
            end
        end
    end
    return res
end

function nexti(i, level)
    (level == 0) && return i + 1
    (level == 1) | (level == 2) | (level == 3) && return i + 8
    (level == 4) | (level == 5) && return i + 4
    (level == 6) | (level == 7) && return i + 2
end

function nextj(j, level)
    (level == 0) | (level == 7) && return j + 1
    (level == 1) | (level == 2) && return j + 8
    (level == 3) | (level == 4) && return j + 4
    (level == 5) | (level == 6) && return j + 2
end

function startpoint(level)
    level == 0 && return (0, 0)   # .+ (1, 1) == (1, 1)
    level == 1 && return (-7, -7) # .+ (8, 8) == (1, 1)
    level == 2 && return (-7, -3) # .+ (8, 8) == (1, 5)
    level == 3 && return (-3, -3) # .+ (8, 4) == (5, 1)
    level == 4 && return (-3, -1) # .+ (4, 4) == (1, 3)
    level == 5 && return (-1, -1) # .+ (4, 2) == (3, 1)
    level == 6 && return (-1, 0)  # .+ (2, 2) == (1, 2)
    level == 7 && return (0, 0)   # .+ (2, 1) == (2, 1)
end

function build!(res::Array{T, 2}, data, idx, width, height, colourtype, bitdepth, plte, level) where T
    i1, j1 = startpoint(level)
    j2 = j1
    @inbounds for i in 1:height
        i1 = nexti(i1, level)
        j1 = j2  # reset counter
        shift = 0
        p = 0x00
        idx += 1
        for j in 1:width
            j1 = nextj(j1, level)
            p2 = if colourtype == 0
                if bitdepth != 16
                    # color bits are written from left to right
                    if shift == 0
                        p = data[idx]
                        idx += 1
                        if bitdepth == 1
                            shift = 7
                        elseif bitdepth == 2
                            shift = 6
                        elseif bitdepth == 4
                            shift = 4
                        elseif bitdepth == 8
                            shift = 0
                        end
                    else
                        if bitdepth == 1
                            shift -= 1
                        elseif bitdepth == 2
                            shift -= 2
                        elseif bitdepth == 4
                            shift -= 4
                        elseif bitdepth == 8
                            shift -= 8
                        end
                    end
                    if bitdepth == 1
                        p2 = (p >> shift) & 0x01 == 0x01 ? 0xff : 0x00
                    elseif bitdepth == 2
                        p2 = ((p >> shift) & 0x03) * 0x55
                    elseif bitdepth == 4
                        p2 = ((p >> shift) & 0x0f) * 0x11
                    elseif bitdepth == 8
                        p2 = p
                    end

                    pixel(T, p2)
                else
                    idx += 2
                    p2 = casttouint16(data, idx - 2)

                    pixel(T, p2)
                end
            elseif colourtype == 2
                # Warning, it is possible that PLTE can be used here
                if bitdepth == 8
                    idx += 3
                    pixel(T, data[idx - 3], data[idx - 2], data[idx - 1])
                elseif bitdepth == 16
                    idx += 6
                    r1 = casttouint16(data, idx - 6)
                    g1 = casttouint16(data, idx - 4)
                    b1 = casttouint16(data, idx - 2)
                    pixel(T, r1, g1, b1)
                end
            elseif colourtype == 3
                if shift == 0
                    p = data[idx]
                    idx += 1
                    if bitdepth == 1
                        shift = 7
                    elseif bitdepth == 2
                        shift = 6
                    elseif bitdepth == 4
                        shift = 4
                    end
                else
                    if bitdepth == 1
                        shift -= 1
                    elseif bitdepth == 2
                        shift -= 2
                    elseif bitdepth == 4
                        shift -= 4
                    end
                end
                if bitdepth == 1
                    p2 = (p >> shift) & 0x01
                elseif bitdepth == 2
                    p2 = (p >> shift) & 0x03
                elseif bitdepth == 4
                    p2 = (p >> shift) & 0x0f
                else
                    p2 = p
                end

                plte[p2 + 1]
            elseif colourtype == 4
                if bitdepth == 8
                    idx += 2
                    pixel(T, data[idx - 2], data[idx - 1])
                else
                    idx += 4
                    c = casttouint16(data, idx - 4)
                    alpha = casttouint16(data, idx - 2)
                    pixel(T, c, alpha)
                end
            else
                if bitdepth == 8
                    idx += 4
                    r1 = data[idx - 4]
                    g1 = data[idx - 3]
                    b1 = data[idx - 2]
                    alpha1 = data[idx - 1]
                    pixel(T, r1, g1, b1, alpha1)
                else
                    idx += 8
                    r1 = casttouint16(data, idx - 8)
                    g1 = casttouint16(data, idx - 6)
                    b1 = casttouint16(data, idx - 4)
                    alpha1 = casttouint16(data, idx - 2)
                    pixel(T, r1, g1, b1, alpha1)
                end
            end

            res[i1, j1] = p2
        end
    end

    return idx
end

end # module


text:
using Abaco
using Test

Abaco.loginit()

# My time zone: Saturday, October 30, 2021 3:43:35 PM GMT+02:00 DST
window = 900
ts = 1635600600
ne = "Civetta"

abaco = abaco_init(interval=window) do ts, ne, name, value, inputs
    @info "[identity] age [$ts]: scope: [$ne] $name = $value"
    @test value == 100.0
end

formula(abaco, "r = x == 0 ? y : 2")

ingest(abaco, ts, ne, "x", 0)

@test_throws Abaco.EvalError ingest(abaco, ts, ne, "y", 100)


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(xgroup=["A", "A", "B", "B"],
     x=["X", "Y", "X", "Y"],
     y=[0.236033, 0.346517, 0.312707, 0.00790928],
     Geom.subplot_grid(Geom.bar))

plot(x=[0.488613, 0.210968, 0.951916],
     y=[0.999905, 0.251662, 0.986666],
     xgroup=["A", "A", "B"],
     ygroup=["X", "Y", "Y"],
     Geom.subplot_grid(Geom.bar))


text:
using YaoCompiler
using Test

@testset "runtime" begin
    include("runtime/locations.jl")
end

@testset "compiler" begin
    include("compiler/parse.jl")
    include("compiler/utils.jl")
    include("compiler/circuit.jl")
end


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

module AzureClusterlessHPC

# If key is set, requires re-compilation
if haskey(ENV, "AZ_BATCH_TASK_WORKING_DIR")
    include("runtime/azureclusterlesshpc_light.jl")
else
    using PyCall, Serialization, JSON, Random, SyntaxTree, Logging
    import Base.fetch, Base.setindex!

    export batch_show, batch_clear, Options, fileinclude
    export delete_pool, delete_container, delete_all_jobs

    # Initiliaze PyCall constants
    const batch = PyNULL()
    const azurequeue = PyNULL()
    const azureblob = PyNULL()
    const serviceprinciple = PyNULL()
    const datetime = PyNULL()
    const batchmodels = PyNULL()
    const azureclusterlesshpc = PyNULL()

    # Variables for collecting expressions
    global __expressions__ = nothing    # collect generic expressions tagged with @batchdef
    global __packages__ = nothing   # collect "using" statements tagged with @batchdef

    # Module initialization
    include("core/parse_credentials.jl")

    function __init__()

        # Add path to Python module to sys env
        sys = pyimport("sys")
        sys.path = cat(sys.path, joinpath(dirname(pathof(AzureClusterlessHPC)), "pyinterface"); dims=1)

        # Load pymodules
        copy!(batch, pyimport_conda("azure.batch", "azure-batch==11.0.0", "conda-forge"))
        copy!(azurequeue, pyimport_conda("azure.storage.queue", "azure-storage-queue==1.4.0", "conda-forge"))
        copy!(azureblob, pyimport_conda("azure.storage.blob", "azure-storage-blob==1.3.1", "conda-forge"))
        copy!(serviceprinciple, pyimport_conda("azure.common.credentials", "azure-common", "conda-forge"))
        copy!(datetime, pyimport("datetime"))
        copy!(batchmodels, batch.models)
        copy!(azureclusterlesshpc, pyimport("azureclusterlesshpc"))

        # Get credentials if environment variable is set
        try
            if haskey(ENV, "CREDENTIALS")
                creds = JSON.parsefile(ENV["CREDENTIALS"])
                typeof(creds) != Vector{Any} && (creds = [creds])
                global __get_credentials__ = creds
            else
                global __get_credentials__ = get_credentials(joinpath(dirname(pathof(AzureClusterlessHPC))[1:end-4], "user_data"))
            end
        catch
            @warn "No credential file specified via ENV[\"CREDENTIALS\"]. All calls to cloud SDKs that require authentication will return nothing."
            global __get_credentials__ = nothing
        end
        global __credentials__ = __get_credentials__

        # Initialize global resource and blobfuture variables
        global __num_accounts__ = ~isnothing(__credentials__) ? length(__credentials__) : 0
        global __resources__ = [[] for i=1:__num_accounts__]

        # Read parameters
        global __params__ = try
            if haskey(ENV, "PARAMETERS")
                user_params = JSON.parsefile(ENV["PARAMETERS"])
            else
                user_params = Dict()
            end
            create_parameter_dict(user_params, default_parameters)
        catch
            nothing
        end

        # Create storage container
        global __container__ = try __params__["_BLOB_CONTAINER"] catch; nothing end # Azure container for temp files
        global __verbose__ = parse(Bool, __params__["_VERBOSE"])

        # Global list of pools (start with no pools)
        global __active_pools__ = Array{Dict}(undef, 0)

        # Batch and blob clients
        global __clients__ = create_clients(__credentials__, batch=true, blob=true)
    end

    # Includes
    include("pyinterface/pyinterface.jl")
    include("pyinterface/pyinterface_test.jl")
    include("core/azureclusterlesshpc_base.jl")
    include("core/futures.jl")
    include("core/batch_controller.jl")
    include("core/batch_macros.jl")
    include("core/batch_shortcuts.jl")
    include("core/batch_fetch.jl")
end
end


text:
using VAWTAC
using Test
using HDF5
# using PyPlot
# close("all")

const path = splitdir(@__FILE__)[1]
# include("$path/../src/VAWTAC.jl")

atol = 1e-6

#Common variables
ntheta = 36

r = 3.0
twist = 0.0
delta = 0.0
af = VAWTAC.readaerodyn("$path/airfoils/NACA_0012_mod.dat")
B = 3
solidity = 0.25
chord = solidity*r/B

Vinf = 1.0
rho = 1.225
mu = 1.7894e-5

Omega = 0.0

# baseline
turbines = Array{VAWTAC.Turbine}(undef,1)
turbines[1] = VAWTAC.Turbine(r, chord, twist, delta, B, af, Omega, 0.0, 0.0)
env = VAWTAC.Environment(Vinf, rho, mu)

@testset "VAWTAC.jl one turbine" begin

    # --------- one turbine -------------

    n = 20
    tsrvec = LinRange(1, 7, n)
    CPvec = zeros(n)
    CTvec = zeros(n)
    Rpvec = zeros(n)
    Tpvec = zeros(n)
    Zpvec = zeros(n)
    thetavec = zeros(n,ntheta)
    for i in eachindex(tsrvec)
        tsr = tsrvec[i]
        turbines[1].Omega = Vinf*tsr/r
        CT, CP, Rp, Tp, Zp, theta = VAWTAC.actuatorcylinder(turbines, env, ntheta)
        CPvec[i] = CP[1]
        CTvec[i] = CT[1]
        Rpvec[i] = Rp[1]
        Tpvec[i] = Tp[1]
        Zpvec[i] = Zp[1]
        thetavec[i,:] = theta
    end

    #write to file
    filename = "$path/single_unit_test_data.h5"
    # HDF5.h5open(filename, "w") do file
    #     HDF5.write(file,"CPvec_old",CPvec)
    #     HDF5.write(file,"CTvec_old",CTvec)
    #     HDF5.write(file,"Rpvec_old",Rpvec)
    #     HDF5.write(file,"Tpvec_old",Tpvec)
    #     HDF5.write(file,"Zpvec_old",Zpvec)
    #     HDF5.write(file,"thetavec_old",thetavec)
    # end

    #read from file
    CPvec_old = HDF5.h5read(filename,"CPvec_old")
    CTvec_old = HDF5.h5read(filename,"CTvec_old")
    Rpvec_old = HDF5.h5read(filename,"Rpvec_old")
    Tpvec_old = HDF5.h5read(filename,"Tpvec_old")
    Zpvec_old = HDF5.h5read(filename,"Zpvec_old")
    thetavec_old = HDF5.h5read(filename,"thetavec_old")

    #test
    @test isapprox(CPvec,CPvec_old;atol)
    @test isapprox(CTvec,CTvec_old;atol)
    @test isapprox(Rpvec,Rpvec_old;atol)
    @test isapprox(Tpvec,Tpvec_old;atol)
    @test isapprox(Zpvec,Zpvec_old;atol)
    @test isapprox(thetavec,thetavec_old;atol)

    # figure()
    # plot(tsrvec, CPvec)
    # xlabel("\$\\lambda\$")
    # ylabel("\$C_p\$")
end

@testset "VAWTAC.jl two turbines" begin
    # -----------  two turbines -----------------
    tsr = 3.5
    Omega = Vinf*tsr/r

    turbines = Array{VAWTAC.Turbine}(undef,2)
    turbines[1] = VAWTAC.Turbine(r, chord, twist, delta, B, af, Omega, 0.0, 0.0)
    turbines[2] = VAWTAC.Turbine(r, chord, twist, delta, B, af, -Omega, 0.0, 2*r)
    CT, CP, Rp, Tp, Zp, theta = VAWTAC.actuatorcylinder(turbines, env, ntheta)

    #write to file
    filename = "$path/dual_unit_test_data.h5"
    # HDF5.h5open(filename, "w") do file
    #     HDF5.write(file,"CP_old",CP)
    #     HDF5.write(file,"CT_old",CT)
    #     HDF5.write(file,"Rp_old",Rp)
    #     HDF5.write(file,"Tp_old",Tp)
    #     HDF5.write(file,"Zp_old",Zp)
    #     HDF5.write(file,"theta_old",theta)
    # end

    #read from file
    CP_old = HDF5.h5read(filename,"CP_old")
    CT_old = HDF5.h5read(filename,"CT_old")
    Rp_old = HDF5.h5read(filename,"Rp_old")
    Tp_old = HDF5.h5read(filename,"Tp_old")
    Zp_old = HDF5.h5read(filename,"Zp_old")
    theta_old = HDF5.h5read(filename,"theta_old")

    #test
    @test isapprox(CP,CP_old;atol)
    @test isapprox(CT,CT_old;atol)
    @test isapprox(Rp,Rp_old;atol)
    @test isapprox(Tp,Tp_old;atol)
    @test isapprox(Zp,Zp_old;atol)
    @test isapprox(theta,theta_old;atol)

    # figure()
    # plot(theta, r*Tp[:, 1])  # r *T = Q
    # plot(theta, r*Tp[:, 2])
    # xlabel("\$\\theta\$")
    # ylabel("Q (torque)")
    # xlim([0, 2*pi])

end

max_stars_repo_name: kevmoor/vawt-ac

max_issues_repo_name: kevmoor/vawt-ac

max_forks_repo_name: kevmoor/vawt-ac


text:
function rnn(Wxh, Whh, Wy, hprev, x, y)
    h = tanh.(Whh * hprev + Wxh * x)
    yhat = Why * h
    cost = sum((yhat .- y) .^ 2.0)
    return h, cost
end

let
    Wxh = randn(4096, 4096); Whh = randn(4096, 4096); Why = randn(128, 4096);
    hprev = randn(4096); x = randn(4096); y = randn(128);
    inputs = [:Wxh=>Wxh, :Whh=>Whh, :Why=>Why, :hprev => hprev, :x => x, :y=>y];
    # at least shouldn't throw an error
    df = xdiff(rnn; ctx=Dict(:cost => :cost), inputs...)
    vals = [val for (name, val) in inputs]
    dvals = df(vals...)
    @test size(dvals[1][1]) == 4096
end


text:
module MarginalsTest

using Test
using ForneyLab

using ForneyLab: generateId, addNode!, associate!, Product

# Integration helper
mutable struct MockNode <: FactorNode
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Int,Interface}

    function MockNode(vars::Vector{Variable}; id=generateId(MockNode))
        n_interfaces = length(vars)
        self = new(id, Vector{Interface}(undef, n_interfaces), Dict{Int,Interface}())
        addNode!(currentGraph(), self)

        for idx = 1:n_interfaces
            self.i[idx] = self.interfaces[idx] = associate!(Interface(self), vars[idx])
        end

        return self
    end
end

@testset "MarginalEntry" begin
    FactorGraph()
    x = Variable()
    nd = MockNode([x])

    # Marginal table should be constructable by hand
    marginal_table = [MarginalEntry(x, [nd.i[1]], Nothing)]
    @test isa(marginal_table, MarginalTable)
end

@testset "marginalTable" begin
    # Construct
    #
    #    a    b
    # [1]--[2]--

    FactorGraph()

    a = Variable()
    b = Variable()

    n1 = MockNode([a])
    n2 = MockNode([a, b])

    marginal_table = marginalTable([a, b])

    @test length(marginal_table) == 2
    @test marginal_table[1].target == a
    @test marginal_table[1].interfaces[1] == n1.i[1]
    @test marginal_table[1].interfaces[2] == n2.i[1]
    @test marginal_table[1].marginal_update_rule == Product

    @test marginal_table[2].target == b
    @test marginal_table[2].interfaces[1] == n2.i[2]
    @test marginal_table[2].marginal_update_rule == Nothing
end

end # module


text:
"""
    _bisection(f, left::T, right::T; max_iters=500, tol=1e-14) -> Union{Nothing, T}

Finds the root `x, f(x) ≈ 0` in interval `[left, right]`, returns nothing if convergence fails.
Errors if the interval is non-bracketing (sign of `f(left)` identical to sign of `f(right)`).
"""
function _bisection(f, left, right; max_iters=500, tol=1e-14)
    for _ in 1:max_iters
        f_left, f_right = f(left), f(right)
        signbit(f_left) == signbit(f_right) && error("Interval became non-bracketing.")

        # Terminate if interval length ~ eps()
        mid = (left + right) / 2
        if left == mid || right == mid
            return mid
        end
        # Terminate if within tol of 0; otherwise, bisect
        f_mid = f(mid)
        if abs(f_mid) < tol
            return mid
        end
        if signbit(f_mid) == signbit(f_left)
            left = mid
        elseif signbit(f_mid) == signbit(f_right)
            right = mid
        end
    end
    return nothing
end


text:
module AssetMgmt

using TimeData
using DataFrames
using DataArrays
using Gadfly
using Dates

## using Datetime
#using NLopt

##############################################################################
##
## Exported methods and types
##
##############################################################################

export # functions and types
UniverseModel,
MuSigmaModel,
DistributionModel,
UniverseEstimate,
MuSigmaUniverse,
InitialStrategy,
MuSigmaStrategy,
MaxMu,
MinSigma,
InvestmentStrategy,
SeparateTurnover
## Investments,
## chkEqualsOne,
## invRet,
## makeWeights

## define basic investment type and type constraints
include("universeModels.jl")
include("universeEstimate.jl")
include("utils.jl")
include("initialStrategies.jl")
include("muSigmaFunctions.jl")

include("constraints.jl")
include("investments.jl")
include("plotting.jl")

## define filters
include("tOv_filters.jl")

## ## define strategies
include("strategies.jl")
## include("markowitz.jl")
## include("moments.jl")

## ## performance analysis
include("performance.jl")
## include("diversification.jl")
## include("grouping.jl")
## include("plotting.jl")

## ## basic utilities
## include("utils.jl")
## include("io.jl")

## ## estimation
## include("estimation.jl")
## include("mu_sigma_estimation.jl")

## ## EGARCH
## include("egarch.jl")

end # module


text:
@testset "CheckVcov" begin
    hrs = exampledata("hrs")
    nt = (data=hrs, esample=trues(size(hrs,1)), aux=trues(size(hrs,1)), vce=Vcov.robust())
    @test checkvcov!(nt...) == NamedTuple()
    nt = merge(nt, (vce=Vcov.cluster(:hhidpn),))
    @test checkvcov!(nt...) == (esample=trues(size(hrs,1)),)

    df = DataFrame(hrs)
    allowmissing!(df, :hhidpn)
    nt = merge(nt, (data=df,))
    @test checkvcov!(nt...) == (esample=trues(size(hrs,1)),)

    @test CheckVcov()((data=hrs, esample=trues(size(hrs,1)), aux=trues(size(hrs,1)))) ==
        (data=hrs, esample=trues(size(hrs,1)), aux=trues(size(hrs,1)))
end

@testset "ParseFEterms" begin
    hrs = exampledata(:hrs)

    ret = parsefeterms!(TermSet())
    @test (ret...,) == (TermSet(), Set{FETerm}(), false)
    ts = TermSet((term(1), term(:male)))
    ret = parsefeterms!(ts)
    @test (ret...,) == (ts, Set{FETerm}(), false)

    ts = TermSet(term(1)+term(:male)+fe(:hhidpn))
    ret = parsefeterms!(ts)
    @test (ret...,) == (TermSet(InterceptTerm{false}(), term(:male)),
        Set(([:hhidpn]=>Symbol[],)), true)

    ts = TermSet(fe(:wave)+fe(:hhidpn))
    ret = parsefeterms!(ts)
    @test (ret...,) == (TermSet(InterceptTerm{false}()),
        Set(([:hhidpn]=>Symbol[], [:wave]=>Symbol[])), true)

    # Verify that no change is made on intercept
    ts = TermSet((term(:male), fe(:hhidpn)&term(:wave)))
    ret = parsefeterms!(ts)
    @test (ret...,) == (TermSet(term(:male)), Set(([:hhidpn]=>[:wave],)), false)

    ts = TermSet((term(:male), fe(:hhidpn)&term(:wave)))
    nt = (xterms=ts,)
    @test ParseFEterms()(nt) == (xterms=TermSet(term(:male)),
        feterms=Set(([:hhidpn]=>[:wave],)), has_fe_intercept=false)
end

@testset "GroupFEterms" begin
    nt = (feterms=Set(([:hhidpn]=>[:wave],)),)
    @test groupfeterms(nt...) === nt
    _byid(GroupFEterms()) == false
    @test GroupFEterms()(nt) === nt
end

@testset "MakeFEs" begin
    hrs = exampledata("hrs")
    @test makefes(hrs, FETerm[]) == (allfes=Dict{FETerm,FixedEffect}(),)

    feterms = [[:hhidpn]=>Symbol[], [:hhidpn,:wave]=>[:male]]
    ret = makefes(hrs, feterms)
    @test ret == (allfes=Dict{FETerm,FixedEffect}(feterms[1]=>FixedEffect(hrs.hhidpn),
        feterms[2]=>FixedEffect(hrs.hhidpn, hrs.wave, interaction=_multiply(hrs, [:male]))),)

    allntargs = NamedTuple[(feterms=Set{FETerm}(),), (feterms=Set(feterms),),
        (feterms=Set(([:hhidpn]=>Symbol[],)),)]
    @test Set(combinedargs(MakeFEs(), allntargs)...) ==
        Set(push!(feterms, [:hhidpn]=>Symbol[]))

    nt = (data=hrs, feterms=feterms)
    @test MakeFEs()(nt) == merge(nt, ret)
end

@testset "CheckFEs" begin
    hrs = exampledata("hrs")
    N = size(hrs, 1)
    nt = (feterms=Set{FETerm}(), allfes=Dict{FETerm,FixedEffect}(),
        esample=trues(N), drop_singletons=true)
    @test checkfes!(nt...) == (esample=nt.esample, fes=FixedEffect[],
        fenames=String[], nsingle=0)

    feterm = [:hhidpn]=>Symbol[]
    allfes = makefes(hrs, [[:hhidpn]=>Symbol[], [:wave]=>[:male]]).allfes
    nt = merge(nt, (feterms=Set((feterm,)), allfes=allfes))
    @test checkfes!(nt...) == (esample=trues(N), fes=[FixedEffect(hrs.hhidpn)],
        fenames=["fe_hhidpn"], nsingle=0)

    # fes are sorted by name
    feterms = Set(([:wave]=>Symbol[], [:hhidpn]=>Symbol[]))
    allfes = makefes(hrs, [[:hhidpn]=>Symbol[], [:wave]=>Symbol[]]).allfes
    nt = merge(nt, (feterms=feterms, allfes=allfes))
    @test checkfes!(nt...) == (esample=trues(N),
        fes=[FixedEffect(hrs.hhidpn), FixedEffect(hrs.wave)],
        fenames=["fe_hhidpn", "fe_wave"], nsingle=0)

    df = DataFrame(hrs)
    df = df[(df.wave.==7).|((df.wave.==8).&(df.wave_hosp.==8)), :]
    N = size(df, 1)
    feterm = [:hhidpn]=>Symbol[]
    allfes = makefes(df, [[:hhidpn]=>Symbol[]]).allfes
    nt = merge(nt, (feterms=Set((feterm,)), allfes=allfes, esample=trues(N)))
    kept = df.wave_hosp.==8
    @test checkfes!(nt...) == (esample=kept, fes=[FixedEffect(df.hhidpn)[kept]],
        fenames=["fe_hhidpn"], nsingle=N-sum(kept))

    df = df[df.wave.==7, :]
    N = size(df, 1)
    allfes = makefes(df, [[:hhidpn]=>Symbol[]]).allfes
    nt = merge(nt, (allfes=allfes, esample=trues(N)))
    @test_throws ErrorException checkfes!(nt...)

    nt = merge(nt, (esample=trues(N),))
    @test_throws ErrorException CheckFEs()(nt)
    nt = merge(nt, (esample=trues(N), drop_singletons=false))
    @test CheckFEs()(nt) == merge(nt, (fes=[FixedEffect(df.hhidpn)],
        fenames=["fe_hhidpn"], nsingle=0))
end

@testset "MakeFESolver" begin
    hrs = exampledata("hrs")
    N = size(hrs, 1)
    fes = FixedEffect[FixedEffect(hrs.hhidpn)]
    fenames = [:fe_hhidpn]
    nt = (fes=fes, weights=uweights(N), default(MakeFESolver())...)
    ret = makefesolver(nt...)
    @test ret.feM isa FixedEffects.FixedEffectSolverCPU{Float64}
    nt = merge(nt, (fes=FixedEffect[],))
    @test makefesolver(nt...) == (feM=nothing,)
    @test MakeFESolver()(nt) == merge(nt, (feM=nothing,))
end

@testset "MakeYXCols" begin
    hrs = exampledata("hrs")
    N = size(hrs, 1)
    t1 = InterceptTerm{true}()
    nt = (data=hrs, weights=uweights(N), esample=trues(N), feM=nothing,
        has_fe_intercept=false, default(MakeYXCols())...)
    ret = makeyxcols(nt..., TermSet(term(:oop_spend), term(1)))
    @test ret.yxterms[term(:oop_spend)] isa ContinuousTerm
    @test ret.yxcols == Dict(ret.yxterms[term(:oop_spend)]=>hrs.oop_spend,
        t1=>ones(N))

    @test ret.nfeiterations === nothing
    @test ret.feconverged === nothing

    # Verify that an intercept will be added if not having one
    ret1 = makeyxcols(nt..., TermSet(term(:oop_spend)))
    @test ret1 == ret
    ret1 = makeyxcols(nt..., TermSet(term(:oop_spend), term(0)))
    @test collect(keys(ret1.yxterms)) == [term(:oop_spend)]
    @test ret1.yxcols == ret.yxcols

    wt = Weights(hrs.rwthh)
    nt = merge(nt, (weights=wt,))
    ret = makeyxcols(nt..., TermSet(term(:oop_spend), term(:riearnsemp), term(:male)))
    @test ret.yxcols == Dict(ret.yxterms[term(:oop_spend)]=>hrs.oop_spend.*sqrt.(wt),
        ret.yxterms[term(:riearnsemp)]=>hrs.riearnsemp.*sqrt.(wt),
        ret.yxterms[term(:male)]=>reshape(hrs.male.*sqrt.(wt), N),
        t1=>reshape(sqrt.(wt), N))

    df = DataFrame(hrs)
    df.riearnsemp[1] = NaN
    nt = merge(nt, (data=df,))
    @test_throws ErrorException makeyxcols(nt..., TermSet(term(:riearnsemp), term(1)))
    df.spouse = convert(Vector{Float64}, df.spouse)
    df.spouse[1] = Inf
    @test_throws ErrorException makeyxcols(nt..., TermSet(term(:oop_spend), term(:spouse)))

    df = DataFrame(hrs)
    x = randn(N)
    df.x = x
    wt = uweights(N)
    fes = [FixedEffect(df.hhidpn)]
    feM = AbstractFixedEffectSolver{Float64}(fes, wt, Val{:cpu}, Threads.nthreads())
    nt = merge(nt, (data=df, weights=wt, feM=feM, has_fe_intercept=true))
    ret = makeyxcols(nt..., TermSet(term(:oop_spend), InterceptTerm{false}(), term(:x)))
    resids = reshape(copy(df.oop_spend), N, 1)
    _feresiduals!(resids, feM, 1e-8, 10000)
    resids .*= sqrt.(wt)
    @test ret.yxcols[ret.yxterms[term(:oop_spend)]] == reshape(resids, N)
    # Verify input data are not modified
    @test df.oop_spend == hrs.oop_spend
    @test df.x == x

    df = DataFrame(hrs)
    esample = df.rwthh.> 0
    N = sum(esample)
    wt = Weights(hrs.rwthh[esample])
    fes = [FixedEffect(df.hhidpn[esample])]
    feM = AbstractFixedEffectSolver{Float64}(fes, wt, Val{:cpu}, Threads.nthreads())
    nt = merge(nt, (data=df, esample=esample, weights=wt, feM=feM, has_fe_intercept=true))
    ret = makeyxcols(nt..., TermSet(term(:oop_spend), InterceptTerm{false}()))
    resids = reshape(df.oop_spend[esample], N, 1)
    _feresiduals!(resids, feM, 1e-8, 10000)
    resids .*= sqrt.(wt)
    @test ret.yxcols == Dict(ret.yxterms[term(:oop_spend)]=>reshape(resids, N))
    @test ret.nfeiterations isa Int
    @test ret.feconverged
    # Verify input data are not modified
    @test df.oop_spend == hrs.oop_spend

    allntargs = NamedTuple[(yterm=term(:oop_spend), xterms=TermSet())]
    @test combinedargs(MakeYXCols(), allntargs) == (TermSet(term(:oop_spend)),)
    push!(allntargs, allntargs[1])
    @test combinedargs(MakeYXCols(), allntargs) == (TermSet(term(:oop_spend)),)
    push!(allntargs, (yterm=term(:riearnsemp),
        xterms=TermSet(InterceptTerm{false}())))
    @test combinedargs(MakeYXCols(), allntargs) ==
        (TermSet(term(:oop_spend), term(:riearnsemp), InterceptTerm{false}()),)
    push!(allntargs, (yterm=term(:riearnsemp), xterms=TermSet(term(:male))))
    @test combinedargs(MakeYXCols(), allntargs) ==
        (TermSet(term(:oop_spend), term(:riearnsemp), InterceptTerm{false}(), term(:male)),)

    nt = merge(nt, (data=df, yterm=term(:oop_spend), xterms=TermSet(InterceptTerm{false}())))
    @test MakeYXCols()(nt) == merge(nt, ret)
end

@testset "MakeTreatCols" begin
    hrs = exampledata("hrs")
    N = size(hrs, 1)
    tr = dynamic(:wave, -1)
    pr = nevertreated(11)
    nt = (data=hrs, treatname=:wave_hosp, treatintterms=TermSet(), feM=nothing,
        weights=uweights(N), esample=trues(N), default(MakeTreatCols())...)
    ret = maketreatcols(nt..., tr.time, Dict(-1=>1), IdDict{ValidTimeType,Int}(11=>1))
    @test size(ret.cells) == (20, 2)
    @test length(ret.rows) == 20
    @test size(ret.treatcells) == (12, 2)
    @test length(ret.treatcols) == length(ret.treatrows) == 12
    @test length(ret.treatweights) == length(ret.treatcounts) == length(ret.treatrows)
    @test ret.cells.wave_hosp == repeat(8:11, inner=5)
    @test ret.cells.wave == repeat(7:11, 4)
    @test ret.rows[end] == findall((hrs.wave_hosp.==11).&(hrs.wave.==11))
    @test ret.treatcells.wave_hosp == repeat(8:10, inner=4)
    rel = ret.cells.wave .- ret.cells.wave_hosp
    @test ret.treatcells.rel == rel[(rel.!=-1).&(ret.cells.wave_hosp.!=11)]
    @test ret.treatrows[1] == findall((hrs.wave_hosp.==8).&(hrs.wave.==8))
    col = convert(Vector{Float64}, (hrs.wave_hosp.==8).&(hrs.wave.==8))
    @test ret.treatcols[1] == col
    @test ret.treatweights == ret.treatcounts
    w = ret.treatweights
    @test all(w[ret.treatcells.wave_hosp.==8].==252)
    @test all(w[ret.treatcells.wave_hosp.==9].==176)
    @test all(w[ret.treatcells.wave_hosp.==10].==163)

    nt = merge(nt, (treatintterms=TermSet(term(:male)),))
    ret1 = maketreatcols(nt..., tr.time, Dict(-1=>1), IdDict{ValidTimeType,Int}(11=>1))
    @test size(ret1.cells) == (40, 3)
    @test length(ret1.rows) == 40
    @test size(ret1.treatcells) == (24, 3)
    @test length(ret1.treatcols) == length(ret1.treatrows) == 24
    @test length(ret1.treatweights) == length(ret1.treatcounts) == length(ret1.treatrows)
    @test ret1.cells.wave_hosp == repeat(8:11, inner=10)
    @test ret1.cells.wave == repeat(repeat(7:11, inner=2), 4)
    @test ret1.cells.male == repeat(0:1, 20)
    @test ret1.rows[end] == findall((hrs.wave_hosp.==11).&(hrs.wave.==11).&(hrs.male.==1))
    @test ret1.treatcells.wave_hosp == repeat(8:10, inner=8)
    rel = ret1.cells.wave .- ret1.cells.wave_hosp
    @test ret1.treatcells.rel == rel[(rel.!=-1).&(ret1.cells.wave_hosp.!=11)]
    @test ret1.treatrows[1] == findall((hrs.wave_hosp.==8).&(hrs.wave.==8).&(hrs.male.==0))
    col1 = convert(Vector{Float64}, (hrs.wave_hosp.==8).&(hrs.wave.==8).&(hrs.male.==0))
    @test ret1.treatcols[1] == col1
    @test ret1.treatweights == ret1.treatcounts

    nt = merge(nt, (cohortinteracted=false, treatintterms=TermSet()))
    ret2 = maketreatcols(nt..., tr.time, Dict(-1=>1), IdDict{ValidTimeType,Int}(11=>1))
    @test ret2.cells[1] == ret.cells[1]
    @test ret2.cells[2] == ret.cells[2]
    @test ret2.rows == ret.rows
    @test size(ret2.treatcells) == (6, 1)
    @test length(ret2.treatcols) == length(ret2.treatrows) == 6
    @test length(ret2.treatweights) == length(ret2.treatcounts) == length(ret2.treatrows)
    @test ret2.treatcells.rel == [-3, -2, 0, 1, 2, 3]
    @test ret2.treatrows[1] == findall((hrs.wave.-hrs.wave_hosp.==-3).&(hrs.wave_hosp.!=11))
    col2 = convert(Vector{Float64}, (hrs.wave.-hrs.wave_hosp.==-3).&(hrs.wave_hosp.!=11))
    @test ret2.treatcols[1] == col2
    @test ret2.treatweights == ret2.treatcounts

    nt = merge(nt, (treatintterms=TermSet(term(:male)),))
    ret3 = maketreatcols(nt..., tr.time, Dict(-1=>1), IdDict{ValidTimeType,Int}(11=>1))
    @test ret3.cells[1] == ret1.cells[1]
    @test ret3.cells[2] == ret1.cells[2]
    @test ret3.cells[3] == ret1.cells[3]
    @test ret3.rows == ret1.rows
    @test size(ret3.treatcells) == (12, 2)
    @test length(ret3.treatcols) == length(ret3.treatrows) == 12
    @test length(ret3.treatweights) == length(ret3.treatcounts) == length(ret3.treatrows)
    @test ret3.treatcells.rel == repeat([-3, -2, 0, 1, 2, 3], inner=2)
    @test ret3.treatcells.male == repeat(0:1, 6)
    @test ret3.treatrows[1] ==
        findall((hrs.wave.-hrs.wave_hosp.==-3).&(hrs.wave_hosp.!=11).&(hrs.male.==0))

    df = DataFrame(hrs)
    esample = df.rwthh.> 0
    N = sum(esample)
    wt = Weights(hrs.rwthh[esample])
    fes = [FixedEffect(df.hhidpn[esample])]
    feM = AbstractFixedEffectSolver{Float64}(fes, wt, Val{:cpu}, Threads.nthreads())
    nt = merge(nt, (data=df, feM=feM, weights=wt, esample=esample,
        treatintterms=TermSet(), cohortinteracted=true))
    ret = maketreatcols(nt..., tr.time, Dict(-1=>1), IdDict{ValidTimeType,Int}(11=>1))
    col = reshape(col[esample], N, 1)
    defaults = (default(MakeTreatCols())...,)
    _feresiduals!(col, feM, defaults[2:3]...)
    @test ret.treatcols[1] == (col.*sqrt.(wt))[:]
    @test ret.treatcounts == [252, 252, 252, 251, 176, 176, 176, 175, 162, 160, 163, 162]
    @test ret.treatweights[1] == 1776173

    allntargs = NamedTuple[(tr=tr, pr=pr)]
    @test combinedargs(MakeTreatCols(), allntargs) ==
        (Dict(-1=>1), IdDict{ValidTimeType,Int}(11=>1))
    push!(allntargs, allntargs[1])
    @test combinedargs(MakeTreatCols(), allntargs) ==
        (Dict(-1=>2), IdDict{ValidTimeType,Int}(11=>2))
    push!(allntargs, (tr=dynamic(:wave, [-1,-2]), pr=nevertreated(10:11)))
    @test combinedargs(MakeTreatCols(), allntargs) ==
        (Dict(-1=>3), IdDict{ValidTimeType,Int}(11=>3))
    push!(allntargs, (tr=dynamic(:wave, [-3]), pr=nevertreated(10)))
    @test combinedargs(MakeTreatCols(), allntargs) ==
        (Dict{Int,Int}(), IdDict{ValidTimeType,Int}())

    allntargs = NamedTuple[(tr=tr, pr=pr), (tr=tr, pr=unspecifiedpr())]
    @test combinedargs(MakeTreatCols(), allntargs) ==
        (Dict(-1=>2), IdDict{ValidTimeType,Int}())

    df.wave = settime(Date.(hrs.wave), Year(1))
    df.wave_hosp = settime(Date.(hrs.wave_hosp), Year(1), start=Date(7))
    ret1 = maketreatcols(nt..., tr.time,
        Dict(-1=>1), IdDict{ValidTimeType,Int}(Date(11)=>1))
    @test ret1.cells[1] == Date.(ret.cells[1])
    @test ret1.cells[2] == Date.(ret.cells[2])
    @test ret1.rows == ret.rows
    @test ret1.treatcells[1] == Date.(ret.treatcells[1])
    @test ret1.treatcells[2] == ret.treatcells[2]
    @test ret1.treatrows == ret.treatrows
    @test ret1.treatcols == ret.treatcols
    @test ret1.treatweights == ret.treatweights
    @test ret1.treatcounts == ret.treatcounts

    rot = ifelse.(isodd.(hrs.hhidpn), 1, 2)
    df.wave = RotatingTimeArray(rot, hrs.wave)
    df.wave_hosp = RotatingTimeArray(rot, hrs.wave_hosp)
    e = rotatingtime((1,1,2), (10,11,11))
    ret2 = maketreatcols(nt..., tr.time,
        Dict(-1=>1), IdDict{ValidTimeType,Int}(c=>1 for c in e))
    @test ret2.cells[1] == sort!(append!((rotatingtime(r, ret.cells[1]) for r in (1,2))...))
    rt = append!((rotatingtime(r, 7:11) for r in (1,2))...)
    @test ret2.cells[2] == repeat(rt, 4)
    @test size(ret2.treatcells) == (20, 2)
    @test sort!(unique(ret2.treatcells[1])) ==
        sort!(append!(rotatingtime(1, 8:9), rotatingtime(2, 8:10)))

    df.wave = settime(Date.(hrs.wave), Year(1), rotation=rot)
    df.wave_hosp = settime(Date.(hrs.wave_hosp), Year(1), start=Date(7), rotation=rot)
    e = rotatingtime((1,1,2), Date.((10,11,11)))
    ret3 = maketreatcols(nt..., tr.time,
        Dict(-1=>1), IdDict{ValidTimeType,Int}(c=>1 for c in e))
    @test ret3.cells[1].time == Date.(ret2.cells[1].time)
    @test ret3.cells[2].time == Date.(ret2.cells[2].time)
    @test ret3.rows == ret2.rows
    @test ret3.treatcells[1].time == Date.(ret2.treatcells[1].time)
    @test ret3.treatcells[2] == ret2.treatcells[2]
    @test ret3.treatrows == ret2.treatrows
    @test ret3.treatcols == ret2.treatcols
    @test ret3.treatweights == ret2.treatweights
    @test ret3.treatcounts == ret2.treatcounts

    nt = merge(nt, (data=hrs, tr=tr, pr=pr))
    @test MakeTreatCols()(nt) == merge(nt, (cells=ret.cells, rows=ret.rows,
        treatcells=ret.treatcells, treatrows=ret.treatrows, treatcols=ret.treatcols,
        treatweights=ret.treatweights, treatcounts=ret.treatcounts))
end

@testset "SolveLeastSquares" begin
    hrs = exampledata("hrs")
    N = size(hrs, 1)
    df = DataFrame(hrs)
    df.t2 = fill(2.0, N)
    t1 = InterceptTerm{true}()
    t0 = InterceptTerm{false}()
    tr = dynamic(:wave, -1)
    pr = nevertreated(11)
    yxterms = Dict([x=>apply_schema(x, schema(x, df), StatisticalModel)
        for x in (term(:oop_spend), t1, term(:t2), t0, term(:male), term(:spouse))])
    yxcols0 = Dict(yxterms[term(:oop_spend)]=>hrs.oop_spend, t1=>ones(N, 1),
        yxterms[term(:male)]=>hrs.male, yxterms[term(:spouse)]=>hrs.spouse)
    col0 = convert(Vector{Float64}, (hrs.wave_hosp.==10).&(hrs.wave.==10))
    col1 = convert(Vector{Float64}, (hrs.wave_hosp.==10).&(hrs.wave.==11))
    col2 = convert(Vector{Float64}, (hrs.wave_hosp.==11).&(hrs.wave.==11))
    treatcells0 = VecColumnTable((wave_hosp=[10, 10, 11], rel=[0, 1, 0]))
    treatcols0 = [col0, col1, col2]
    treatcounts0 = [sum(c.==1) for c in treatcols0]
    treatweights0 = convert(Vector{Float64}, treatcounts0)
    nt = (tr=tr, pr=pr, yterm=term(:oop_spend), xterms=TermSet(term(1)), yxterms=yxterms,
        yxcols=yxcols0, treatcells=treatcells0, treatcols=treatcols0,
        treatweights=treatweights0, treatcounts=treatcounts0,
        cohortinteracted=true, has_fe_intercept=false)
    ret = solveleastsquares!(nt...)
    # Compare estimates with Stata
    # gen col0 = wave_hosp==10 & wave==10
    # gen col1 = wave_hosp==10 & wave==11
    # reg oop_spend col0 col1
    @test ret.coef[1] ≈ 2862.4141 atol=1e-4
    @test ret.coef[2] ≈ 490.44869 atol=1e-4
    @test ret.coef[3] ≈ 3353.6565 atol=1e-4
    @test ret.treatcells.wave_hosp == [10, 10]
    @test ret.treatcells.rel == [0, 1]
    @test ret.xterms == AbstractTerm[t1]
    @test ret.basiscols == trues(3)

    # Verify that an intercept will only be added when needed
    nt1 = merge(nt, (xterms=TermSet(),))
    ret1 = solveleastsquares!(nt1...)
    @test ret1 == ret
    nt1 = merge(nt, (xterms=TermSet(term(0)),))
    ret1 = solveleastsquares!(nt1...)
    @test ret1.xterms == AbstractTerm[]
    @test size(ret1.X, 2) == 2
    nt1 = merge(nt, (xterms=TermSet((term(:spouse), term(:male))),))
    ret1 = solveleastsquares!(nt1...)
    @test ret1.xterms == AbstractTerm[yxterms[term(:male)], yxterms[term(:spouse)],
        InterceptTerm{true}()]

    # Test colliner xterms are handled
    yxcols1 = Dict(yxterms[term(:oop_spend)]=>hrs.oop_spend, t1=>ones(N),
        yxterms[term(:t2)]=>df.t2)
    nt1 = merge(nt, (xterms=TermSet(term(1), term(:t2)), yxcols=yxcols1))
    ret1 = solveleastsquares!(nt1...)
    @test ret1.coef[1:2] == ret.coef[1:2]
    @test sum(ret1.basiscols) == 3

    treatcells1 = VecColumnTable((rel=[0, 1, 1, 1], wave_hosp=[10, 10, 0, 1]))
    treatcols1 = push!(treatcols0[1:2], ones(N), ones(N))
    treatweights1 = push!(copy(treatweights0), 10.0)
    treatcounts1 = push!(copy(treatcounts0), 10)
    # basecol is rather conservative in dropping collinear columns
    # If there are three constant columns, it may be that only one of them gets dropped
    # Also need to have at least one term in xterms for basiscol to work
    yxcols2 = Dict(yxterms[term(:oop_spend)]=>hrs.oop_spend, yxterms[term(:male)]=>hrs.male)
    nt1 = merge(nt1, (xterms=TermSet(term(:male), term(0)), treatcells=treatcells1,
        yxcols=yxcols2, treatcols=treatcols1, treatweights=treatweights1,
        treatcounts=treatcounts1))
    @test_throws ErrorException solveleastsquares!(nt1...)

    @test SolveLeastSquares()(nt) == merge(nt, ret)
end

@testset "EstVcov" begin
    hrs = exampledata("hrs")
    N = size(hrs, 1)
    col0 = convert(Vector{Float64}, (hrs.wave_hosp.==10).&(hrs.wave.==10))
    col1 = convert(Vector{Float64}, (hrs.wave_hosp.==10).&(hrs.wave.==11))
    y = convert(Vector{Float64}, hrs.oop_spend)
    X = hcat(col0, col1, ones(N, 1))
    crossx = cholesky!(Symmetric(X'X))
    cf = crossx \ (X'y)
    residuals = y - X * cf
    nt = (data=hrs, esample=trues(N), vce=Vcov.simple(), coef=cf, X=X, crossx=crossx,
        residuals=residuals, xterms=AbstractTerm[term(1)], fes=FixedEffect[],
        has_fe_intercept=false)
    ret = estvcov(nt...)
    # Compare estimates with Stata
    # reg oop_spend col0 col1
    # mat list e(V)
    @test ret.vcov_mat[1,1] ≈ 388844.2 atol=0.1
    @test ret.vcov_mat[2,2] ≈ 388844.2 atol=0.1
    @test ret.vcov_mat[3,3] ≈ 20334.169 atol=1e-3
    @test ret.vcov_mat[2,1] ≈ 20334.169 atol=1e-3
    @test ret.vcov_mat[3,1] ≈ -20334.169 atol=1e-3
    @test ret.dof_resid == N - 3
    @test ret.F ≈ 10.68532285556941 atol=1e-6

    nt = merge(nt, (vce=Vcov.robust(), fes=FixedEffect[]))
    ret = estvcov(nt...)
    # Compare estimates with Stata
    # reg oop_spend col0 col1, r
    # mat list e(V)
    @test ret.vcov_mat[1,1] ≈ 815817.44 atol=0.1
    @test ret.vcov_mat[2,2] ≈ 254993.93 atol=1e-2
    @test ret.vcov_mat[3,3] ≈ 19436.209 atol=1e-3
    @test ret.vcov_mat[2,1] ≈ 19436.209 atol=1e-3
    @test ret.vcov_mat[3,1] ≈ -19436.209 atol=1e-3
    @test ret.dof_resid == N - 3
    @test ret.F ≈ 5.371847047691197 atol=1e-6

    nt = merge(nt, (vce=Vcov.cluster(:hhidpn),))
    ret = estvcov(nt...)
    # Compare estimates with Stata
    # reghdfe oop_spend col0 col1, noa clu(hhidpn)
    # mat list e(V)
    @test ret.vcov_mat[1,1] ≈ 744005.01 atol=0.1
    @test ret.vcov_mat[2,2] ≈ 242011.45 atol=1e-2
    @test ret.vcov_mat[3,3] ≈ 28067.783 atol=1e-3
    @test ret.vcov_mat[2,1] ≈ 94113.386 atol=1e-2
    @test ret.vcov_mat[3,1] ≈ 12640.559 atol=1e-2
    @test ret.dof_resid == N - 3
    @test ret.F ≈ 5.542094561672688 atol=1e-6

    fes = FixedEffect[FixedEffect(hrs.hhidpn)]
    wt = uweights(N)
    feM = AbstractFixedEffectSolver{Float64}(fes, wt, Val{:cpu}, Threads.nthreads())
    X = hcat(col0, col1)
    _feresiduals!(Combination(y, X), feM, 1e-8, 10000)
    crossx = cholesky!(Symmetric(X'X))
    cf = crossx \ (X'y)
    residuals = y - X * cf
    nt = merge(nt, (vce=Vcov.robust(), coef=cf, X=X, crossx=crossx, residuals=residuals,
        xterms=AbstractTerm[], fes=fes, has_fe_intercept=true))
    ret = estvcov(nt...)
    # Compare estimates with Stata
    # reghdfe oop_spend col0 col1, a(hhidpn) vce(robust)
    # mat list e(V)
    @test ret.vcov_mat[1,1] ≈ 654959.97 atol=0.1
    @test ret.vcov_mat[2,2] ≈ 503679.27 atol=0.1
    @test ret.vcov_mat[2,1] ≈ 192866.2 atol=0.1
    @test ret.dof_resid == N - nunique(fes[1]) - 2
    @test ret.F ≈ 7.559815337537517 atol=1e-6

    nt = merge(nt, (vce=Vcov.cluster(:hhidpn),))
    ret = estvcov(nt...)
    # Compare estimates with Stata
    # reghdfe oop_spend col0 col1, a(hhidpn) clu(hhidpn)
    # mat list e(V)
    @test ret.vcov_mat[1,1] ≈ 606384.66 atol=0.1
    @test ret.vcov_mat[2,2] ≈ 404399.89 atol=0.1
    @test ret.vcov_mat[2,1] ≈ 106497.43 atol=0.1
    @test ret.dof_resid == N - 3
    @test ret.F ≈ 8.197452252592386 atol=1e-6

    @test EstVcov()(nt) == merge(nt, ret)
end

@testset "SolveLeastSquaresWeights" begin
    hrs = exampledata("hrs")
    tr = dynamic(:wave, -1)
    N = size(hrs, 1)
    col0 = convert(Vector{Float64}, (hrs.wave_hosp.==10).&(hrs.wave.==10))
    col1 = convert(Vector{Float64}, (hrs.wave_hosp.==10).&(hrs.wave.==11))
    cellnames = [:wave_hosp, :wave]
    cols = subcolumns(hrs, cellnames)
    cells, rows = cellrows(cols, findcell(cols))
    yterm = term(:oop_spend)
    yxterms = Dict(x=>apply_schema(x, schema(x, hrs), StatisticalModel)
        for x in (term(:oop_spend), term(:male)))

    yxcols = Dict(yxterms[term(:oop_spend)]=>convert(Vector{Float64}, hrs.oop_spend))
    wt = uweights(N)
    fes = [FixedEffect(hrs.hhidpn)]
    feM = AbstractFixedEffectSolver{Float64}(fes, wt, Val{:cpu}, Threads.nthreads())
    y = yxcols[yxterms[term(:oop_spend)]]
    fetol = 1e-8
    _feresiduals!(Combination(y, col0, col1), feM, fetol, 10000)
    X = hcat(col0, col1)
    crossx = cholesky!(Symmetric(X'X))
    cf = crossx \ (X'y)
    tcells = VecColumnTable((rel=[0, 1],))

    nt = (tr=tr, solvelsweights=true, lswtnames=(), cells=cells, rows=rows,
        X=X, crossx=crossx, coef=cf, treatcells=tcells, yterm=yterm, xterms=AbstractTerm[],
        yxterms=yxterms, yxcols=yxcols, feM=feM, fetol=1e-8, femaxiter=10000, weights=wt)
    ret = solveleastsquaresweights(nt...)
    lswt = ret.lsweights
    @test lswt.r === cells
    @test lswt.c === tcells
    @test lswt[lswt.r.wave_hosp.==10, 1] ≈ [-1/3, -1/3, -1/3, 1, 0] atol=fetol
    @test lswt[lswt.r.wave_hosp.==10, 1] ≈ [-1/3, -1/3, -1/3, 1, 0] atol=fetol
    @test lswt[lswt.r.wave_hosp.==10, 2] ≈ [-1/3, -1/3, -1/3, 0, 1] atol=fetol
    @test all(x->isapprox(x, 0, atol=fetol), lswt[lswt.r.wave_hosp.!=10, :])
    @test ret.cellweights == ret.cellcounts == length.(rows)
    @test all(i->ret.cellymeans[i] ≈ sum(y[rows[i]])/length(rows[i]), 1:length(rows))

    nt0 = merge(nt, (lswtnames=(:no,),))
    @test_throws ArgumentError solveleastsquaresweights(nt0...)

    nt = merge(nt, (lswtnames=(:wave,),))
    ret = solveleastsquaresweights(nt...)
    lswt = ret.lsweights
    @test size(lswt.r) == (5, 1)
    @test lswt.r.wave == 7:11

    nt = merge(nt, (lswtnames=(:wave, :wave_hosp),))
    ret = solveleastsquaresweights(nt...)
    lswt = ret.lsweights
    @test size(lswt.r) == (20, 2)
    @test lswt.r.wave == repeat(7:11, inner=4)
    @test lswt.r.wave_hosp == repeat(8:11, outer=5)
    @test lswt[lswt.r.wave_hosp.==10, 1] ≈ [-1/3, -1/3, -1/3, 1, 0] atol=fetol
    @test lswt[lswt.r.wave_hosp.==10, 2] ≈ [-1/3, -1/3, -1/3, 0, 1] atol=fetol

    df = DataFrame(hrs)
    df.x = zeros(N)
    df.x[df.wave.==7] .= rand(656)
    yxterms = Dict(x=>apply_schema(x, schema(x, df), StatisticalModel)
        for x in (term(:oop_spend), term(:x)))
    yxcols = Dict(yxterms[term(:oop_spend)]=>convert(Vector{Float64}, hrs.oop_spend),
        yxterms[term(:x)]=>convert(Vector{Float64}, df.x))
    fes = [FixedEffect(hrs.hhidpn), FixedEffect(hrs.wave)]
    feM = AbstractFixedEffectSolver{Float64}(fes, wt, Val{:cpu}, Threads.nthreads())
    y = yxcols[yxterms[term(:oop_spend)]]
    x = yxcols[yxterms[term(:x)]]
    _feresiduals!(Combination(y, col0, col1, x), feM, fetol, 10000)
    X = hcat(col0, col1, x)
    crossx = cholesky!(Symmetric(X'X))
    cf = crossx \ (X'y)
    xterms = AbstractTerm[yxterms[term(:x)]]
    nt = merge(nt, (lswtnames=(:wave_hosp, :wave), X=X, crossx=crossx, coef=cf,
        xterms=xterms, yxterms=yxterms, yxcols=yxcols))
    ret = solveleastsquaresweights(nt...)
    lswt = ret.lsweights
    @test lswt[lswt.r.wave.<=9, 1] ≈ lswt[lswt.r.wave.<=9, 2]
    @test lswt[lswt.r.wave.==10, 1] ≈ lswt[lswt.r.wave.==11, 2]
    @test lswt[lswt.r.wave.==11, 1] ≈ lswt[lswt.r.wave.==10, 2]
    y1 = y .- cf[3].*x
    @test all(i->ret.cellymeans[i] == sum(y1[rows[i]])/length(rows[i]), 1:length(rows))

    @test SolveLeastSquaresWeights()(nt) == merge(nt, ret)
end


text:
"""
    unionof(T₁, T₂, ...)

Similar to `Union{T₁, T₂, ...}` but the order is preserved. This object can
be passed to `UnionArray` to specify the element type.
"""
function UnionArrays.unionof(types::Type...)
    vals = foldlargs((), types...) do vals, T
        (vals..., Val{T}())
    end
    return verify(UnionOf(vals))
end

const NTypes{N} = NTuple{N, Val}

struct UnionOf{T<:NTypes}
    types::T
end

const ElTypeSpec = Union{Type, NTypes, UnionOf}

astupletype(u::UnionOf) =
    foldrargs(u.types..., Tuple{}) do v, T
        Base.tuple_type_cons(valueof(v), T)
    end

asunion(u::UnionOf) = asunion(u.types)

asunion(types::NTypes) =
    foldlargs(Union{}, types...) do T, v
        Union{T,valueof(v)}
    end

asntypes(types::NTypes) = types
asntypes(u::UnionOf) = u.types
@generated asntypes(::Type{T}) where {T} =
    QuoteNode(foldrunion((S, types) -> (Val(S), types...), T, ()))

function verify(u::UnionOf)
    if length(u.types) != union_ntypes(asunion(u))
        error("non-unique types are provided")
    end
    return u
end

function Base.show(io::IO, u::UnionOf)
    @nospecialize u
    print(io, UnionArrays.unionof)
    print(io, '(')
    isfirst = true
    for v in u.types
        if !isfirst
            print(io, ", ")
        end
        isfirst = false
        print(io, valueof(v))
    end
    print(io, ')')
end


text:
"""
A database transaction. Every operation requires a transaction handle.
All database operations require a transaction handle. Transactions may be read-only or read-write.
"""
type Transaction
    handle::Ptr{Void}
    Transaction() = new(C_NULL)
    Transaction(h::Ptr{Void}) = new(h)
end

function env(txn::Transaction)
    env_ptr = ccall((:mdb_txn_env, liblmdb), Ptr{Void}, (Ptr{Void},), txn.handle)
    (env_ptr == C_NULL) && return nothing
    return Environment(env_ptr)
end

"Check if transaction is open."
isopen(txn::Transaction) = txn.handle != C_NULL

"""Create a transaction for use with the environment

`start` function creates a new transaction and returns `Transaction` object.
It allows to set transaction flags with `flags` option.
"""
function start(env::Environment; flags::Cuint=zero(Cuint),
               parent::Nullable{Transaction} = Nullable{Transaction}())
    txn_ref = Ref{Ptr{Void}}(C_NULL)
    ret = ccall( (:mdb_txn_begin, liblmdb), Cint,
                  (Ptr{Void}, Ptr{Void}, Cuint, Ptr{Ptr{Void}}),
                   env.handle, get(parent, Transaction()).handle,  flags, txn_ref)
    (ret != 0) && throw(LMDBError(ret))
    return Transaction(txn_ref[])
end
start(f::Function, env::Environment; flags::Cuint=zero(Cuint)) = f(start(env, flags=flags))

"""Abandon all the operations of the transaction instead of saving them

The transaction and its cursors must not be used after, because its handle is freed.
"""
function abort(txn::Transaction)
    ccall( (:mdb_txn_abort, liblmdb), Void, (Ptr{Void},), txn.handle)
    txn.handle = C_NULL
    return
end

"""Commit all the operations of a transaction into the database

The transaction and its cursors must not be used after, because its handle is freed.
"""
function commit(txn::Transaction)
    ret = ccall( (:mdb_txn_commit, liblmdb), Cint, (Ptr{Void},), txn.handle)
    txn.handle = C_NULL
    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"""Reset a read-only transaction

Abort the transaction like `abort`, but keep the transaction handle.
"""
function reset(txn::Transaction)
    ret = ccall( (:mdb_txn_reset, liblmdb), Cint, (Ptr{Void},), txn.handle)
    (ret != 0) && throw(LMDBError(ret))
    return ret
end

"""Renew a read-only transaction

This acquires a new reader lock for a transaction handle that had been released by `reset`.
It must be called before a reset transaction may be used again.
"""
function renew(txn::Transaction)
    ret = ccall( (:mdb_txn_renew, liblmdb), Cint, (Ptr{Void},), txn.handle)
    (ret != 0) && throw(LMDBError(ret))
    return ret
end


text:
module PP2PM
export load_pm_from_json

import JSON

function load_pm_from_json(json_path)
    pm = Dict()
    open(json_path, "r") do f
        pm = JSON.parse(f)  # parse and transform data
    end

    for (idx, gen) in pm["gen"]
        if gen["model"] == 1
            pm["gen"][idx]["cost"] = convert(Array{Float64,1}, gen["cost"])
        end
    end
    return pm
end

end

max_stars_repo_name: Zamwell/pandapower

max_issues_repo_name: Zamwell/pandapower

max_forks_repo_name: Zamwell/pandapower


text:
using TOML

function get_fasta_info( input_filename )

    line = readlines( input_filename )

    i = 1
    while i ≤ length( line )
        if line[ i ][ begin ] == '>'
            push!( title_list, line[ i ] )
            push!( seq_list, line[ i + 1 ] )
            i += 1
        else
            seq_list[ end ] *= line[ i ]
        end
        i += 1
    end

end

function check_query_exist( query_title, input_filename )

    query_length = -1
    for i = 1 : length( title_list )
        if ( title_list )[ i ]  == query_title
            query_length = length( seq_list[ i ] )
            println( "\nQuery sequence :" )
            println( query_title, " (", query_length, " AA)" )
            break
        end
    end

    if query_length == -1
        println( "ERROR !" )
        println( "The query sequence does not exist in '", input_filename, " ' : " )
        println( query_title )
        println( "Program halted !!!" )
        exit( 1 )
    end

    return query_length

end

function check_threshold( short_threshold, long_threshold )

    if !( 0 ≤ short_threshold ≤ 1 )
        println( "ERROR !" )
        println( "Lower sequence length threshold must be bounded with a range of zero to one." )
        println( "Program halted !!!" )
        exit( 1 )
    end

    if !( long_threshold ≥ 1 )
        println( "ERROR !" )
        println( "Upper sequence length threshold must be more than one." )
        println( "Program halted !!!" )
        exit( 1 )
    end

end

function remove_short_long( short_length, long_length )

    println( "\nRemoved sequences : " )
    println( "Shorter than ", short_length, " AA or longer than ", long_length, " AA.\n" )

    for i = 1 : length( title_list )
        if length( seq_list[ i ] ) < short_length || long_length < length( seq_list[ i ] )
            println( title_list[ i ], " (",  length( seq_list[ i ] ), " AA)" )
        else
            push!( output_seq_list, seq_list[ i ] )
            push!( output_title_list, title_list[ i ] )
        end
    end

end

function save_result( output_filename )

    open( output_filename, "w") do fout

        for i = 1 : length( output_title_list )
            write( fout, output_title_list[ i ], "\n" )
            for j = 1 : length( output_seq_list[ i ] )
                if j % 60 == 0
                    write( fout, output_seq_list[ i ][ j ], "\n" )
                else
                    write( fout, output_seq_list[ i ][ j ] )
                end
            end
            if length( output_seq_list[ i ] ) % 60 != 0
                write( fout, "\n")
            end
        end

    end;

    println( "\n'", output_filename, "' was correctly written.\n" )
end

println( "\nRemove too short and too long sequences compared to the query sequence in MSA.\n" )

# Check the number of arguments.
if length( ARGS ) != 1
    println( "Only one argument, input filename." )
    exit( 1 )
end

seq_list          = []
title_list        = []
output_seq_list   = []
output_title_list = []

# Read TOML file and get input information.
input = read( ARGS[ 1 ], String )
#println( input )
toml_data = TOML.parse( input )

input_filename  = toml_data[ "file" ][ "input" ][ "name" ]
output_filename = toml_data[ "file" ][ "output" ][ "name" ]
query_title     = toml_data[ "query" ][ "title" ]
short_threshold = toml_data[ "threshold" ][ "short" ][ "value" ]
long_threshold  = toml_data[ "threshold" ][ "long" ][ "value" ]

# Read Multi-FASTA file and get sequence information.
get_fasta_info( input_filename )

# Check whether or not the query sequence exists in MSA.
query_length = check_query_exist( query_title, input_filename )

# Check thresholds are correct.
check_threshold( short_threshold, long_threshold )

short_length = Int( floor( query_length * short_threshold ) )
long_length  = Int( floor( query_length * long_threshold ) )

# Remove too short and long sequences.
remove_short_long( short_length, long_length )

# Save result.
save_result( output_filename )

println( "Program completed !!!" )

max_stars_repo_name: shin-kinos/rem-diff-len

max_issues_repo_name: shin-kinos/rem-diff-len

max_forks_repo_name: shin-kinos/rem-diff-len


text:
using JsonDB
using Test

const d = JsonDB

@testset "JsonDB.jl" begin
	db = d.load("test.db", false)

	@test db !== nothing

	d.set(db, "fizz", "buzz")
	@test "buzz" == d.get(db, "fizz")
	@test "buzz" == db.db["fizz"]

	d.deletedb(db)
	d.set(db, "key1", "value1")
	d.set(db, "key2", "value2")
	d.create_dict(db, "dict1")
	d.create_array(db, "array1")
	output = d.getall(db)
	expected = collect(keys(Dict("key2"=>nothing, "key1"=>nothing, "dict1"=>nothing, "array1"=>nothing)))
	@test collect(output) == expected

	d.set(db, "key", "value")
	d.delete(db, "key")
	@test d.get(db, "key") == false

	d.set(db, "key", "value")
	d.append(db, "key", "value")
	@test d.get(db, "key") == "valuevalue"

	d.set(db, "key", "value")
	@test haskey(db, "key") == true
	d.delete(db, "key")
	@test haskey(db, "key") == false

	d.create_array(db, "list")
	d.set_array(db, "list", "one")
	d.set_array(db, "list", "two")
	d.set_array(db, "list", "three")
	d.set_array(db, "list", "four")
	@test d.range(db, "list", 2, 3) == ["two", "three"]
	d.delete_array(db, "list")

	d.create_dict(db, "dict")
	d.set_dict(db, "dict", "key" => "value")
	@test d.haskey(db, "dict", "key") == true
	d.delete_dict(db, "dict")

	d.create_dict(db, "dict")
	d.set_dict(db, "dict", "key" => "value")
	@test d.haskey(db, "dict", "absent_key") == false
	d.delete_dict(db, "dict")
end


text:
#
# Dataspaces
#

const DS_SCALAR = 0x00
const DS_SIMPLE = 0x01
const DS_NULL = 0x02

struct WriteDataspace{N,A<:Tuple}
    dataspace_type::UInt8
    size::NTuple{N,Length}
    attributes::A
end

struct ReadDataspace
    dataspace_type::UInt8
    dimensionality::UInt8
    dimensions_offset::Int64
end
ReadDataspace() = ReadDataspace(DS_SCALAR, 0, -1)

struct DataspaceStart
    version::UInt8
    dimensionality::UInt8
    flags::UInt8
    dataspace_type::UInt8
end
define_packed(DataspaceStart)

const EMPTY_DIMENSIONS = Int[]

# Pass-through for custom serializations
# Need a bunch of methods to avoid ambiguity
WriteDataspace(f::JLDFile, x, ::Type{CustomSerialization{S,ODR}}) where {S,ODR} =
    WriteDataspace(f, x, ODR)
WriteDataspace(f::JLDFile, x::Array, ::Type{CustomSerialization{S,ODR}}) where {S,ODR} =
    WriteDataspace(f, x, ODR)
WriteDataspace(f::JLDFile, x::Array{T,0}, ::Type{CustomSerialization{S,ODR}}) where {T,S,ODR} =
    WriteDataspace(f, x, ODR)

WriteDataspace() = WriteDataspace(DS_NULL, (), ())
WriteDataspace(::JLDFile, ::Any, odr::Nothing) = WriteDataspace()
WriteDataspace(::JLDFile, ::Any, ::Any) = WriteDataspace(DS_SCALAR, (), ())

# Ghost type array
WriteDataspace(f::JLDFile, x::Array{T}, ::Nothing) where {T} =
   WriteDataspace(DS_NULL, (),
             (WrittenAttribute(f, :dimensions, collect(Int64, reverse(size(x)))),))

# Reference array
WriteDataspace(f::JLDFile, x::Array{T,N}, ::Type{RelOffset}) where {T,N} =
    WriteDataspace(DS_SIMPLE, convert(Tuple{Vararg{Length}}, reverse(size(x))),
              (WrittenAttribute(f, :julia_type, write_ref(f, T, f.datatype_wsession)),))

# isbitstype array
WriteDataspace(f::JLDFile, x::Array, ::Any) =
    WriteDataspace(DS_SIMPLE, convert(Tuple{Vararg{Length}}, reverse(size(x))), ())

# Zero-dimensional arrays need an empty dimensions attribute
WriteDataspace(f::JLDFile, x::Array{T,0}, ::Nothing) where {T} =
    WriteDataspace(DS_NULL, (Length(1),),
              (WrittenAttribute(f, :dimensions, EMPTY_DIMENSIONS)))
WriteDataspace(f::JLDFile, x::Array{T,0}, ::Type{RelOffset}) where {T} =
    WriteDataspace(DS_SIMPLE, (Length(1),),
              (WrittenAttribute(f, :julia_type, write_ref(f, T, f.datatype_wsession)),
               WrittenAttribute(f, :dimensions, EMPTY_DIMENSIONS)))
WriteDataspace(f::JLDFile, x::Array{T,0}, ::Any) where {T} =
    WriteDataspace(DS_SIMPLE, (Length(1),),
              (WrittenAttribute(f, :dimensions, EMPTY_DIMENSIONS),))

jlsizeof(::Union{WriteDataspace{N},Type{WriteDataspace{N}}}) where {N} = 4 + jlsizeof(Length)*N
numel(x::WriteDataspace{0}) = x.dataspace_type == DS_SCALAR ? 1 : 0
numel(x::WriteDataspace) = Int(prod(x.size))

function jlwrite(io::IO, dspace::WriteDataspace{N}) where N
    jlwrite(io, DataspaceStart(2, N, 0, dspace.dataspace_type))
    for x in dspace.size
        jlwrite(io, x::Length)
    end
end

# Reads a dataspace from the file and returns a
# (dataspace_type::UInt8, dataspace_dimensions:G:Vector{Length})
# tuple, where dataspace_type is one of the DS_* constants and
# dataspace_dimensions are the corresponding dimensions
function read_dataspace_message(io::IO)
    dspace_start = jlread(io, DataspaceStart)
    dspace_start.version == 2 || throw(UnsupportedVersionException())
    dataspace_type = dspace_start.dataspace_type
    ReadDataspace(dataspace_type, dspace_start.dimensionality, position(io))
end


text:
import Bijectors: OrderedBijector

@testset "OrderedBijector" begin
    b = OrderedBijector()

    # Length 1
    x = randn(1)
    y = b(x)
    test_bijector(b, hcat(x, x), hcat(y, y), zeros(2))

    # Larger
    x = randn(5)
    xs = hcat(x, x)
    test_bijector(b, xs)

    y = b(x)
    @test sort(y) == y

    ys = b(xs)
    @test sort(ys[:, 1]) == ys[:, 1]
    @test sort(ys[:, 2]) == ys[:, 2]
end