text:
module MaxHelpingHandDirectionalJumpThrus

using ..Ahorn, Maple

# ==== Upside-down jumpthrus

@mapdef Entity "MaxHelpingHand/UpsideDownJumpThru" UpsideDownJumpThru(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, texture::String="wood",
    animationDelay::Number=0.0, pushPlayer::Bool=false, attached::Bool=false)

const textures = ["wood", "dream", "temple", "templeB", "cliffside", "reflection", "core", "moon"]

const quads = Tuple{Integer, Integer, Integer, Integer}[
    (0, 0, 8, 7) (8, 0, 8, 7) (16, 0, 8, 7);
    (0, 8, 8, 5) (8, 8, 8, 5) (16, 8, 8, 5)
]

Ahorn.editingOptions(entity::UpsideDownJumpThru) = Dict{String, Any}(
    "texture" => textures
)

Ahorn.minimumSize(entity::UpsideDownJumpThru) = 8, 0
Ahorn.resizable(entity::UpsideDownJumpThru) = true, false

function Ahorn.selection(entity::UpsideDownJumpThru)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::UpsideDownJumpThru, room::Maple.Room)
    texture = get(entity.data, "texture", "wood")
    texture = texture == "default" ? "wood" : texture

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))

    startX = div(x, 8) + 1
    stopX = startX + div(width, 8) - 1
    startY = div(y, 8) + 1
    animated = Number(get(entity.data, "animationDelay", 0)) > 0

    Ahorn.Cairo.save(ctx)

    Ahorn.scale(ctx, 1, -1)

    len = stopX - startX
    for i in 0:len
        if animated
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)00", 8 * i, -8)
        else
            connected = false
            qx = 2
            if i == 0
                connected = get(room.fgTiles.data, (startY, startX - 1), false) != '0'
                qx = 1

            elseif i == len
                connected = get(room.fgTiles.data, (startY, stopX + 1), false) != '0'
                qx = 3
            end

            quad = quads[2 - connected, qx]
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)", 8 * i, -8, quad...)
        end
    end

    Ahorn.Cairo.restore(ctx)
end

function Ahorn.flipped(entity::UpsideDownJumpThru, horizontal::Bool)
    if !horizontal
        return RegularJumpThru(entity.x, entity.y, entity.width, entity.texture)
    end
end

Ahorn.rotated(entity::UpsideDownJumpThru, steps::Int) = SidewaysJumpThru(entity.x, entity.y, entity.width, steps > 0, entity.texture)

# ==== regular jumpthru got from the result of rotating or flipping helping hand jumpthrus

@mapdef Entity "MaxHelpingHand/RegularJumpThru" RegularJumpThru(x::Integer, y::Integer, width::Integer=8, texture::String="wood", surfaceIndex::Int16=convert(Int16, -1))

Ahorn.editingOptions(entity::RegularJumpThru) = Dict{String, Any}(
    "texture" => textures,
    "surfaceIndex" => Maple.tileset_sound_ids
)

Ahorn.minimumSize(entity::RegularJumpThru) = 8, 0
Ahorn.resizable(entity::RegularJumpThru) = true, false

function Ahorn.selection(entity::RegularJumpThru)
    x, y = Ahorn.position(entity)
    width = Int(get(entity.data, "width", 8))

    return Ahorn.Rectangle(x, y, width, 8)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::RegularJumpThru, room::Maple.Room)
    texture = get(entity.data, "texture", "wood")
    texture = texture == "default" ? "wood" : texture

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))

    startX = div(x, 8) + 1
    stopX = startX + div(width, 8) - 1
    startY = div(y, 8) + 1

    len = stopX - startX
    for i in 0:len
        connected = false
        qx = 2
        if i == 0
            connected = get(room.fgTiles.data, (startY, startX - 1), false) != '0'
            qx = 1

        elseif i == len
            connected = get(room.fgTiles.data, (startY, stopX + 1), false) != '0'
            qx = 3
        end

        quad = quads[2 - connected, qx]
        Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)", 8 * i, 0, quad[1], quad[2], quad[3], quad[4])
    end
end

function Ahorn.flipped(entity::RegularJumpThru, horizontal::Bool)
    if !horizontal
        return UpsideDownJumpThru(entity.x, entity.y, entity.width, entity.texture)
    end
end

Ahorn.rotated(entity::RegularJumpThru, steps::Int) = SidewaysJumpThru(entity.x, entity.y, entity.width, steps < 0, entity.texture)

# ==== Sideways jumpthrus

@mapdef Entity "MaxHelpingHand/SidewaysJumpThru" SidewaysJumpThru(x::Integer, y::Integer, height::Integer=Maple.defaultBlockHeight,
    left::Bool=true, texture::String="wood", animationDelay::Number=0.0, letSeekersThrough::Bool=false, surfaceIndex::Int=-1, pushPlayer::Bool=false)
@mapdef Entity "MaxHelpingHand/AttachedSidewaysJumpThru" AttachedSidewaysJumpThru(x::Integer, y::Integer, height::Integer=Maple.defaultBlockHeight,
    left::Bool=true, texture::String="wood", animationDelay::Number=0.0, letSeekersThrough::Bool=false, surfaceIndex::Int=-1, pushPlayer::Bool=false)

const placements = Ahorn.PlacementDict(
    "Upside Down Jump Through (max480's Helping Hand)" => Ahorn.EntityPlacement(
        UpsideDownJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood"
        )
    ),
    "Sideways Jump Through (Left) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        SidewaysJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => true
        )
    ),
    "Sideways Jump Through (Right) (max480's Helping Hand)" => Ahorn.EntityPlacement(
        SidewaysJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => false
        )
    ),
    "Sideways Jump Through (Attached, Left)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        AttachedSidewaysJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => true
        )
    ),
    "Sideways Jump Through (Attached, Right)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        AttachedSidewaysJumpThru,
        "rectangle",
        Dict{String, Any}(
            "texture" => "wood",
            "left" => false
        )
    )
)

# Sideways Jumpthrus and Attached ones look the same and have the same options, so we're going to handle them together!
const jumpthruUnion = Union{SidewaysJumpThru, AttachedSidewaysJumpThru}

Ahorn.editingOptions(entity::jumpthruUnion) = Dict{String, Any}(
    "texture" => textures,
    "surfaceIndex" => Maple.tileset_sound_ids
)

function Ahorn.flipped(entity::jumpthruUnion, horizontal::Bool)
    if horizontal
        entity.left = !entity.left
        return entity
    end
end

function Ahorn.rotated(entity::jumpthruUnion, steps::Int)
    if (steps > 0) == entity.left
        return RegularJumpThru(entity.x, entity.y, entity.height, entity.texture, convert(Int16, entity.surfaceIndex))
    else
        return UpsideDownJumpThru(entity.x, entity.y, entity.height, entity.texture)
    end
end

Ahorn.minimumSize(entity::jumpthruUnion) = 0, 8
Ahorn.resizable(entity::jumpthruUnion) = false, true

function Ahorn.selection(entity::jumpthruUnion)
    x, y = Ahorn.position(entity)
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, 8, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::jumpthruUnion, room::Maple.Room)
    texture = get(entity.data, "texture", "wood")
    texture = texture == "default" ? "wood" : texture

    # Values need to be system specific integer
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    height = Int(get(entity.data, "height", 8))
    left = get(entity.data, "left", true)

    startX = div(x, 8) + 1
    startY = div(y, 8) + 1
    stopY = startY + div(height, 8) - 1
    animated = Number(get(entity.data, "animationDelay", 0)) > 0

    Ahorn.Cairo.save(ctx)

    Ahorn.rotate(ctx, pi / 2)

    if left
        Ahorn.scale(ctx, 1, -1)
    end

    len = stopY - startY
    for i in 0:len
        if animated
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)00", 8 * i, left ? 0 : -8)
        else
            connected = false
            qx = 2
            if i == 0
                connected = get(room.fgTiles.data, (startY - 1, startX), false) != '0'
                qx = 1

            elseif i == len
                connected = get(room.fgTiles.data, (stopY + 1, startX), false) != '0'
                qx = 3
            end

            quad = quads[2 - connected, qx]
            Ahorn.drawImage(ctx, "objects/jumpthru/$(texture)", 8 * i, left ? 0 : -8, quad...)
        end
    end

    Ahorn.Cairo.restore(ctx)
end

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
"""
    Dimension

Abstract supertype of all dimension types.

Example concrete implementations are [`X`](@ref), [`Y`](@ref), [`Z`](@ref),
[`Ti`](@ref) (Time), and the custom [`Dim`]@ref) dimension.

`Dimension`s label the axes of an [`AbstractDimArray`](@ref),
or other dimensional objects, and are used to index into the array.

They may also provide an alternate index to lookup for each array axis.
This may be any `AbstractVector` matching the array axis length, or a `Val`
holding a tuple for compile-time index lookups.

`Dimension`s also have `mode` and `metadata` fields.

`mode` gives more details about the dimension, such as that it is
[`Categorical`](@ref) or [`Sampled`](@ref) as [`Points`](@ref) or
[`Intervals`](@ref) along some transect. DimensionalData will
attempt to guess the mode from the passed-in index value.

`metadata` can hold any metadata object adding more information about
the array axis - useful for extending DimensionalData for specific
contexts, like geospatial data in GeoData.jl. By default it is `NoMetadata()`.

Example:

```jldoctest Dimension
using DimensionalData, Dates

x = X(2:2:10)
y = Y(['a', 'b', 'c'])
ti = Ti(DateTime(2021, 1):Month(1):DateTime(2021, 12))

A = DimArray(zeros(3, 5, 12), (y, x, ti))

# output

3×5×12 DimArray{Float64,3} with dimensions:
  Y: Char[a, b, c] (Categorical - Unordered)
  X: 2:2:10 (Sampled - Ordered Regular Points)
  Ti (Time): DateTime("2021-01-01T00:00:00"):Month(1):DateTime("2021-12-01T00:00:00") (Sampled - Ordered Regular Points)
[:, :, 1]
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
[and 11 more slices...]
```

For simplicity, the same `Dimension` types are also used as wrappers
in `getindex`, like:

```jldoctest Dimension
x = A[X(2), Y(3)]

# output

12-element DimArray{Float64,1} with dimensions:
  Ti (Time): DateTime("2021-01-01T00:00:00"):Month(1):DateTime("2021-12-01T00:00:00") (Sampled - Ordered Regular Points)
and reference dimensions:
  Y: c (Categorical - Unordered)
  X: 4 (Sampled - Ordered Regular Points)
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
```

A `Dimension` can also wrap [`Selector`](@ref).

```jldoctest Dimension
x = A[X(Between(3, 4)), Y(At('b'))]

# output

1×12 DimArray{Float64,2} with dimensions:
  X: 4:2:4 (Sampled - Ordered Regular Points)
  Ti (Time): DateTime("2021-01-01T00:00:00"):Month(1):DateTime("2021-12-01T00:00:00") (Sampled - Ordered Regular Points)
and reference dimensions:
  Y: b (Categorical - Unordered)
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```

`Dimension` objects may have [`mode`](@ref) and [`metadata`](@ref) fields
to track additional information about the data and the index, and their relationship.
"""
abstract type Dimension{T,Mo,Me} end

"""
    IndependentDim <: Dimension

Abstract supertype for independent dimensions. Thise will plot on the X axis.
"""
abstract type IndependentDim{T,Mo,Me} <: Dimension{T,Mo,Me} end

"""
    DependentDim <: Dimension

Abstract supertype for Dependent dimensions. These will plot on the Y axis.
"""
abstract type DependentDim{T,Mo,Me} <: Dimension{T,Mo,Me} end

"""
    XDim <: IndependentDim

Abstract supertype for all X dimensions.
"""
abstract type XDim{T,Mo,Me} <: IndependentDim{T,Mo,Me} end

"""
    YDim <: DependentDim

Abstract supertype for all Y dimensions.
"""
abstract type YDim{T,Mo,Me} <: DependentDim{T,Mo,Me} end

"""
    ZDim <: DependentDim

Abstract supertype for all Z dimensions.
"""
abstract type ZDim{T,Mo,Me} <: DependentDim{T,Mo,Me} end

"""
    TimeDim <: IndependentDim

Abstract supertype for all time dimensions.

In a `TimeDime` with `Interval` sampling the locus will automatically
be set to `Start()`. Dates and times generally refer to the start of a
month, hour, second etc., not the central point as is more common with spatial data.
`"""
abstract type TimeDim{T,Mo,Me} <: IndependentDim{T,Mo,Me} end

ConstructionBase.constructorof(d::Type{<:Dimension}) = basetypeof(d)
Adapt.adapt_structure(to, dim::Dimension) =
    rebuild(dim, adapt(to, val(dim)), adapt(to, mode(dim)), adapt(to, metadata(dim)))

const DimType = Type{<:Dimension}
const DimTuple = Tuple{<:Dimension,Vararg{<:Dimension}}
const DimTypeTuple = Tuple{<:DimType,Vararg{<:DimType}}
const VectorOfDim = Vector{<:Dimension}
const DimOrDimType = Union{Dimension,DimType}
const AllDims = Union{Dimension,DimTuple,DimType,DimTypeTuple,VectorOfDim}

# DimensionalData interface methods

"""
    rebuild(dim::Dimension, val, mode=mode(dim), metadata=metadata(dim)) => Dimension
    rebuild(dim::Dimension; val=val(dim), mode=mode(dim), metadata=metadata(dim)) => Dimension

Rebuild dim with fields from `dim`, and new fields passed in.
"""
function rebuild(
    dim::D, val, mode::IndexMode=mode(dim), metadata=metadata(dim)
) where D <: Dimension
    constructorof(D)(val, mode, metadata)
end

dims(dim::Union{Dimension,DimType,Val{<:Dimension}}) = dim
dims(dims::DimTuple) = dims
dims(x) = nothing
dims(::Nothing) = error("No dims found")

refdims(x) = ()

val(dim::Dimension) = dim.val
mode(dim::Dimension) = dim.mode
mode(dim::Union{DimType,Val{<:Dimension}}) = NoIndex()
metadata(dim::Dimension) = dim.metadata

index(dim::Dimension{<:AbstractArray}) = val(dim)
index(dim::Dimension{<:Val}) = unwrap(val(dim))

name(dim::Dimension) = name(typeof(dim))
name(dim::Val{D}) where D = name(D)

bounds(dim::Dimension) = _bounds(mode(dim), dim)

modetype(dim::Dimension) = typeof(mode(dim))
modetype(::Type{<:Dimension{<:Any,Mo}}) where Mo = Mo
modetype(x) = NoIndex

for func in (:order, :span, :sampling, :locus)
    @eval ($func)(dim::Dimension) = ($func)(mode(dim))
end

# Dipatch on Tuple{<:Dimension}, and map to single dim methods
for f in (:val, :index, :mode, :metadata, :order, :sampling, :span, :bounds, :locus,
             :name, :label, :units)
    @eval begin
        $f(dims::DimTuple) = map($f, dims)
        $f(dims::Tuple{}) = ()
        $f(dims::DimTuple, l1, lookup...) = $f(dims, (l1, lookup...))
        $f(dims::DimTuple, lookup) = $f(DD.dims(dims, key2dim(lookup)))
    end
end

order(ot::Type{<:SubOrder}, dims::DimTuple) = map(d -> order(ot, d), dims)
order(ot::Type{<:SubOrder}, dims::Tuple{}) = ()
order(ot::Type{<:SubOrder}, dims_::DimTuple, l1, l2, ls...) =
    order(ot, dims_, (l1, l2, ls...))
order(ot::Type{<:SubOrder}, dims_::DimTuple, lookup::Tuple) =
    map(d -> order(ot, d), dims(dims_, key2dim(lookup)))
order(ot::Type{<:SubOrder}, dims_::DimTuple, lookup) =
    order(ot, dims(dims_, key2dim(lookup)))
order(ot::Type{<:SubOrder}, dim::Dimension) = order(ot, mode(dim))

# Base methods
const ArrayOrVal = Union{AbstractArray,Val}

Base.eltype(d::Type{<:Dimension{T}}) where T = T
Base.eltype(d::Type{<:Dimension{A}}) where A<:AbstractArray{T} where T = T
# TODO Use a vector wrapper instead of Val
Base.eltype(d::Type{<:Dimension{<:Val{Index}}}) where Index = typeof(first(Index))
Base.size(d::Dimension) = size(val(d))
Base.size(d::Dimension{<:Val}) = (length(index(d)),)
Base.axes(d::Dimension{<:AbstractArray}) = axes(index(d))
Base.axes(d::Dimension{<:Val}) = (Base.OneTo(length(d)),)
Base.axes(d::Dimension, i) = axes(index(d), i)
Base.eachindex(d::Dimension{<:Number}) = eachindex(val(d))
Base.eachindex(d::Dimension{<:ArrayOrVal}) = eachindex(index(d))
Base.length(d::Dimension{<:Number}) = length(val(d))
Base.length(d::Dimension{<:ArrayOrVal}) = length(index(d))
Base.ndims(d::Dimension) = 0
Base.ndims(d::Dimension{<:AbstractArray}) = ndims(val(d))
Base.ndims(d::Dimension{<:Val}) = 1
@inline Base.getindex(d::Dimension{<:Number}) = val(d)
@propagate_inbounds Base.getindex(d::Dimension{<:AbstractArray}, i) =
    getindex(index(d), sel2indices(d, i))
@propagate_inbounds Base.getindex(d::Dimension{<:Val{Index}}, i) where Index =
    getindex(Index, sel2indices(d, i))

Base.iterate(d::Dimension{<:ArrayOrVal}, args...) = iterate(index(d), args...)
Base.first(d::Dimension) = val(d)
Base.first(d::Dimension{<:ArrayOrVal}) = first(index(d))
Base.last(d::Dimension) = val(d)
Base.last(d::Dimension{<:ArrayOrVal}) = last(index(d))
Base.firstindex(d::Dimension) = 1
Base.lastindex(d::Dimension) = 1
Base.firstindex(d::Dimension{<:ArrayOrVal}) = firstindex(index(d))
Base.lastindex(d::Dimension{<:ArrayOrVal}) = lastindex(index(d))
Base.step(d::Dimension) = step(mode(d), d)
Base.Array(d::Dimension{<:ArrayOrVal}) = collect(index(d))
function Base.:(==)(d1::Dimension, d2::Dimension)
    typeof(d1) == typeof(d2) &&
    val(d1) == val(d2) &&
    mode(d1) == mode(d2) &&
    metadata(d1) == metadata(d2)
end

"""
Abstract supertype for Dimensions with user-set type paremeters
"""
abstract type ParametricDimension{X,T,Mo,Me} <: Dimension{T,Mo,Me} end

"""
    Dim{S}()
    Dim{S}(val=:; mode=AutoMode(), metadata=NoMetadata())
    Dim{S}(val, mode, metadata=NoMetadata())

A generic dimension. For use when custom dims are required when loading
data from a file. Can be used as keyword arguments for indexing.

Dimension types take precedence over same named `Dim` types when indexing
with symbols, or e.g. creating Tables.jl keys.

```jldoctest
using DimensionalData

dim = Dim{:custom}(['a', 'b', 'c'])

# output

Dim{:custom}:
  val: Char[a, b, c]
Dim{:custom, Vector{Char}, AutoMode{AutoOrder}, NoMetadata}
```
"""
struct Dim{S,T,Mo<:Mode,Me<:AllMetadata} <: ParametricDimension{S,T,Mo,Me}
    val::T
    mode::Mo
    metadata::Me
    Dim{S}(val::T, mode::Mo, metadata::Me=NoMetadata()) where {S,T,Mo<:Mode,Me} =
        new{S,T,Mo,Me}(val, mode, metadata)
end
Dim{S}(val=:; mode::Mode=AutoMode(), metadata=NoMetadata()) where S =
    Dim{S}(val, mode, metadata)

name(::Type{<:Dim{S}}) where S = S
basetypeof(::Type{<:Dim{S}}) where S = Dim{S}
key2dim(s::Val{S}) where S = Dim{S}()
dim2key(::Type{D}) where D<:Dim{S} where S = S

"""
    AnonDim <: Dimension

    AnonDim()

Anonymous dimension. Used when extra dimensions are created,
such as during transpose of a vector.
"""
struct AnonDim{T} <: Dimension{T,NoIndex,NoMetadata}
    val::T
end
AnonDim() = AnonDim(Colon())
AnonDim(val, arg1, args...) = AnonDim(val)

mode(::AnonDim) = NoIndex()
metadata(::AnonDim) = NoMetadata()
name(::AnonDim) = :Anon

"""
    @dim typ [supertype=Dimension] [name::String=string(typ)]

Macro to easily define new dimensions. The supertype will be inserted
into the type of the dim. The default is simply `YourDim <: Dimension`. Making
a Dimesion inherit from `XDim`, `YDim`, `ZDim` or `TimeDim` will affect
automatic plot layout and other methods that dispatch on these types. `<: YDim`
are plotted on the Y axis, `<: XDim` on the X axis, etc.

Example:
```julia
@dim Lat YDim "latitude"
@dim Lon XDim "Longitude"
```
"""
macro dim end
macro dim(typ::Symbol, args...)
    dimmacro(typ::Symbol, :(DimensionalData.Dimension), args...)
end
macro dim(typ::Symbol, supertyp::Symbol, args...)
    dimmacro(typ, supertyp, args...)
end

function dimmacro(typ, supertype, name::String=string(typ))
    quote
        Base.@__doc__ struct $typ{T,Mo<:DimensionalData.Mode,Me<:DimensionalData.AllMetadata} <: $supertype{T,Mo,Me}
            val::T
            mode::Mo
            metadata::Me
        end
        function $typ(
             val=:; mode::DimensionalData.Mode=DimensionalData.AutoMode(),
             metadata=DimensionalData.NoMetadata()
        )
            $typ(val, mode, metadata)
        end
        function $typ(val::V, mode::Mo) where {V,Mo<:DimensionalData.Mode}
            $typ{V,Mo,NoMetadata}(val, mode, DimensionalData.NoMetadata())
        end
        DimensionalData.name(::Type{<:$typ}) = $(QuoteNode(Symbol(name)))
        DimensionalData.key2dim(::Val{$(QuoteNode(typ))}) = $typ()
    end |> esc
end

# Define some common dimensions.

"""
    X <: XDim

    X(val=:; mode=AutoMode(), metadata=NoMetadata())

X [`Dimension`](@ref). `X <: XDim <: IndependentDim`

## Example:
```julia
xdim = X(2:2:10)
# Or
val = A[X(1)]
# Or
mean(A; dims=X)
```
"""
@dim X XDim

"""
    Y <: YDim

    Y(val=:; mode=AutoMode(), metadata=NoMetadata())

Y [`Dimension`](@ref). `Y <: YDim <: DependentDim`

## Example:
```julia
ydim = Y(['a', 'b', 'c'])
# Or
val = A[Y(1)]
# Or
mean(A; dims=Y)
```
"""
@dim Y YDim

"""
    Z <: ZDim

    Z(val=:; mode=AutoMode(), metadata=NoMetadata())

Z [`Dimension`](@ref). `Z <: ZDim <: Dimension`

## Example:
```julia
zdim = Z(10:10:100)
# Or
val = A[Z(1)]
# Or
mean(A; dims=Z)
```
"""
@dim Z ZDim

"""m
    Ti <: TimeDim

    Ti(val=:; mode=AutoMode(), metadata=NoMetadata())

Time [`Dimension`](@ref). `Ti <: TimeDim <: IndependentDim`

`Time` is already used by Dates, and `T` is a common type parameter,
We use `Ti` to avoid clashes.

## Example:
```julia
timedim = Ti(DateTime(2021, 1):Month(1):DateTime(2021, 12))
# Or
val = A[Ti(1)]
# Or
mean(A; dims=Ti)
```
"""
@dim Ti TimeDim "Time"

# Time dimensions need to default to the Start() locus, as that is
# nearly always the format and Center intervals are difficult to
# calculate with DateTime step values.
identify(locus::AutoLocus, dimtype::Type{<:TimeDim}, index) = Start()

const Time = Ti # For some backwards compat

"""
    formatdims(A, dims) => Tuple{Vararg{<:Dimension,N}}

Format the passed-in dimension(s) `dims` to match the array `A`.

This means converting indexes of `Tuple` to `LinRange`, and running
`identify`. Errors are also thrown if dims don't match the array dims or size.

If a [`IndexMode`](@ref) hasn't been specified, an mode is chosen
based on the type and element type of the index:
"""
formatdims(A::AbstractArray, dims) = formatdims(A, (dims,))
formatdims(A::AbstractArray, dims::NamedTuple) = begin
    dims = map((k, v) -> Dim{k}(v), keys(dims), values(dims))
    formatdims(axes(A), dims)
end
formatdims(A::AbstractArray{<:Any,N}, dims::Tuple{Vararg{<:Any,N}}) where N =
    formatdims(axes(A), dims)
@noinline formatdims(A::AbstractArray{<:Any,N}, dims::Tuple{Vararg{<:Any,M}}) where {N,M} =
    throw(DimensionMismatch("Array A has $N axes, while the number of dims is $M"))
formatdims(axes::Tuple, dims::Tuple) = _formatdims(axes, dims)

_formatdims(axes::Tuple{Vararg{<:AbstractRange}}, dims::Tuple) =
    map(_formatdims, axes, dims)
_formatdims(axis::AbstractRange, dimname::Symbol) = Dim{dimname}(axis, NoIndex(), NoMetadata())
_formatdims(axis::AbstractRange, T::Type{<:Dimension}) = T(axis, NoIndex(), NoMetadata())
_formatdims(axis::AbstractRange, dim::Dimension) = begin
    checkaxis(dim, axis)
    rebuild(dim, val(dim), identify(mode(dim), basetypeof(dim), val(dim)))
end
_formatdims(axis::AbstractRange, dim::Dimension{<:NTuple{2}}) = begin
    start, stop = val(dim)
    range = LinRange(start, stop, length(axis))
    _formatdims(axis, rebuild(dim, range))
end
# Dimensions holding colon dispatch on mode
_formatdims(axis::AbstractRange, dim::Dimension{Colon}) = _formatdims(mode(dim), axis, dim)
_formatdims(mode::Union{AutoMode,NoIndex}, axis::AbstractRange, dim::Dimension{Colon}) =
    rebuild(dim, axis, NoIndex())
_formatdims(mode::IndexMode, axis::AbstractRange, dim::Dimension{Colon}) =
    rebuild(dim, axis, identify(mode, basetypeof(dim), axis))

function checkaxis(dim, axis)
    if !(first(axes(dim)) == axis)
        throw(DimensionMismatch(
            "axes of $(basetypeof(dim)) of $(first(axes(dim))) do not match array axis of $axis"
        ))
    end
end

ext: jl


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#

module Raytracing

using Base: Bool, String, Int64
using Colors, Images, ImageIO, FileIO, Polynomials, Test, ArgParse
using ColorTypes:RGB
using LinearAlgebra, StaticArrays
using Printf, ProgressMeter
using Documenter, DocStringExtensions, JSON, Dates,  Intervals

import Base.:+; import Base.:-; import Base.:≈; import Base.:/; import Base.:*
import Base: write, read, print, println, copy, length;

import LinearAlgebra.:⋅; import LinearAlgebra.:×

# from Structs.jl
export BLACK, WHITE, mutable_for_test, to_RGB, HDRimage, get_matrix, Parameters
export Point, Vec, Normal, VEC_X, VEC_Y, VEC_Z, Transformation
export Ray, Camera, OrthogonalCamera, PerspectiveCamera, ImageTracer
export Pigment, UniformPigment, CheckeredPigment, ImagePigment
export BRDF, DiffuseBRDF, SpecularBRDF, Material
export Shape, Sphere, Plane, Torus, Triangle, Cube, AABB
export Vec2d, HitRecord, PointLight, World
export Renderer, OnOffRenderer, FlatRenderer, PathTracer, PointLightRenderer
# from Operations.jl
export are_close, squared_norm, norm, normalize #, normalized_dot
# from ReadingWriting.jl
export load_image, ldr2pfm
# from RangeTesters.jl
export check_is_positive, check_is_uint64, check_is_even_uint64
export check_is_square, check_is_color, check_is_vector
export check_is_declare_float, check_is_one_of, check_is_iterable
export check_is_function, check_is_vec_variables
# from ParseSettings.jl
export parse_command_line, parse_demo_settings
export parse_tonemapping_settings, parse_demoanimation_settings
export parse_render_settings, parse_render_animation_settings
# from ToneMapping.jl
export normalize_image!,  clamp_image!, γ_correction!, tone_mapping
# from Transformations.jl
export rotation_x, rotation_y, rotation_z
export scaling, translation, inverse, is_consistent
# from ImageTracer.jl
export at, fire_ray, fire_all_rays!
# from Shapes.jl
export ray_intersection, add_shape!, add_light!
export is_point_visible, quick_ray_intersection
# from Demo.jl
export demo, demo_animation
# from Pigment.jl
export get_color, evaluate
# from PCG.jl
export PCG, random
# from OrthoNormalBasis.jl
export create_onb_from_z
#from ScatterRay.jl
export scatter_ray
# from Interpreter.jl
# see the module Interpreter.jl
# from Render.jl
export render
# from RenderAnimation.jl
export render_animation
# from PrintFunctions.jl
export print_not_black
# from ArgParse_CLI.jl
export ArgParse_command_line, print_ArgParseSettings

CAMERAS = ["ort", "per"]
RENDERERS = ["onoff", "flat", "pathtracer", "pointlight"]
DEMO_WORLD_TYPES = ["A", "B"]

SYM_NUM = Dict("e"=>ℯ, "pi"=>π)

export CAMERAS, RENDERERS, DEMO_WORLD_TYPES

include("PCG.jl")
include("Structs.jl")
include("Operations.jl")
include("ReadingWriting.jl")
include("RangeTesters.jl")
include("ParseSettings.jl")
include("ToneMapping.jl")
include("Transformations.jl")
include("ImageTracer.jl")
include("Shapes.jl")
include("Demo.jl")
include("Pigment.jl")
include("Renderers.jl")
include("OrthoNormalBasis.jl")
include("ScatterRay.jl")

include("Interpreter.jl")
using .Interpreter

include("Render.jl")
include("RenderAnimation.jl")
include("PrintFunctions.jl")
include("ArgParse_CLI.jl")

# Here is were you can define your own functions to be used in
# the scene-files!
include("YOUR_FUNCTIONS.jl")

end


text:
import Base: read, write, readavailable
import DandelionWebSockets:
    start, start_reader, stop, WriterTaskProxy, FrameFromServer, SocketClosed, ClientProtocolTaskProxy

network_test_frame4 =
    Frame(true, false, false, false, OPCODE_BINARY, false, 126, 256, nomask, zero256)

facts("Reader task") do
    context("Start and stop") do
        # Start async reader task
        # Check that it's running.
        s = BlockingStream(IOBuffer())
        logic = MockClientProtocol([
            (:handle, Any[SocketClosed()])
        ])

        logic_proxy = ClientProtocolTaskProxy(logic)
        start(logic_proxy)

        @sync begin
            @async begin
                reader = start_reader(s, logic_proxy)
                sleep(0.3)
                @fact reader.task --> istaskstarted
                @fact reader.task --> not(istaskdone)

                # Stop reader task
                # Check that it isn't running.
                stop(reader)
                sleep(0.3)
                @fact istaskdone(reader.task) --> true
            end
        end

        check(logic)
    end

    context("Start and stop, with one frame read") do
        framebuf = IOBuffer(Array{UInt8, 1}(), true, true)
        write(framebuf, test_frame1)
        iobuf = IOBuffer(take!(framebuf))
        s = BlockingStream(iobuf)
        logic = MockClientProtocol([
            (:handle, Any[FrameFromServer(test_frame1)]),
            (:handle, Any[SocketClosed()])
        ])

        logic_proxy = ClientProtocolTaskProxy(logic)
        start(logic_proxy)

        @sync begin
            @async begin
                reader = start_reader(s, logic_proxy)
                sleep(0.3)
                @fact reader.task --> istaskstarted
                @fact reader.task --> not(istaskdone)

                # Stop reader task
                # Check that it isn't running.
                stop(reader)
                sleep(0.3)
                @fact istaskdone(reader.task) --> true
            end
        end

        check(logic)
    end

    context("Read a frame and stop") do
        framebuf = IOBuffer(Array{UInt8, 1}(), true, true)
        write(framebuf, test_frame1)
        iobuf = IOBuffer(take!(framebuf))
        s = BlockingStream(iobuf)

        logic = MockClientProtocol([
            (:handle, Any[FrameFromServer(test_frame1)]),
            (:handle, Any[SocketClosed()])
        ])

        logic_proxy = ClientProtocolTaskProxy(logic)
        start(logic_proxy)

        @sync begin
            @async begin
                reader = start_reader(s, logic_proxy)
                sleep(0.3)
                stop(reader)
                sleep(0.3)
                @fact reader.task --> istaskdone
            end
        end
    end

    context("Read several frames and stop") do
        framebuf = IOBuffer(Array{UInt8, 1}(), true, true)
        write(framebuf, test_frame1)
        write(framebuf, test_frame2)
        write(framebuf, test_frame3)
        iobuf = IOBuffer(take!(framebuf))
        s = BlockingStream(iobuf)

        logic = MockClientProtocol([
            (:handle, Any[FrameFromServer(test_frame1)]),
            (:handle, Any[FrameFromServer(test_frame2)]),
            (:handle, Any[FrameFromServer(test_frame3)]),
            (:handle, Any[SocketClosed()])
        ])

        logic_proxy = ClientProtocolTaskProxy(logic)
        start(logic_proxy)

        @sync begin
            @async begin
                reader = start_reader(s, logic_proxy)
                sleep(0.005)
                stop(reader)
                sleep(0.005)
                @fact reader.task --> istaskdone
            end
        end
    end
end

type MockFrameStream <: IO
    frames::Vector{Frame}
end

Base.write(s::MockFrameStream, frame::Frame) = push!(s.frames, frame)

function expect(s::MockFrameStream, frame::Frame)
    @fact s.frames --> x -> !isempty(x)

    actual = shift!(s.frames)
    @fact actual --> frame
end

facts("Writer task") do
    context("Stop and start writer") do
        s = MockFrameStream(Vector{Frame}())
        writer = WriterTaskProxy(s)
        task = start(writer)
        sleep(0.05)
        @fact task --> not(istaskdone)

        stop(writer)
        sleep(0.05)
        @fact task --> istaskdone
    end

    context("Write a few frames") do
        s = MockFrameStream(Vector{Frame}())

        writer = WriterTaskProxy(s)
        task = start(writer)
        sleep(0.05)
        @fact task --> not(istaskdone)

        write(writer, test_frame1)
        write(writer, test_frame2)
        write(writer, test_frame3)

        sleep(0.1)

        expect(s, test_frame1)
        expect(s, test_frame2)
        expect(s, test_frame3)

        stop(writer)
        sleep(0.1)
        @fact task --> istaskdone
    end
end


text:
# differentiated version of sbp_sat_shock_capturing.jl

function calcShockCapturing_diff(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts,
                             capture::SBPParabolicReduced2SC{Tsol, Tres},
                             shockmesh::ShockedElements,
                             assem::AssembleElementData) where {Tsol, Tres}

  println("\nEntered calcShockCapturing_diff")
  computeGradW(mesh, sbp, eqn, opts, capture, shockmesh,
               capture.entropy_vars, capture.diffusion)

  computeVolumeTerm_diff(mesh, sbp, eqn, opts, capture,
                               capture.diffusion, capture.entropy_vars,
                               shockmesh, assem)

  computeFaceTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
                       capture.diffusion, capture.entropy_vars,
                       capture.penalty, assem)

#  computeDirichletBoundaryTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
#                                    assem)

#  if shockmesh.isNeumann
#    println("computing Neumann boundary condition")
#    computeNeumannBoundaryTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
#                              capture.diffusion, capture.entropy_vars, assem)
#  else
#    println("computing Dirichlet boundary condition")
#    computeDirichletBoundaryTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
#                                      assem)
#  end

  #@time computeSharedFaceTerm_diff(mesh, sbp, eqn, opts, capture, shockmesh,
  #                                 capture.diffusion, capture.penalty)

  return nothing
end

function computeVolumeTerm_diff(mesh, sbp, eqn, opts,
                           capture::SBPParabolicReduced2SC{Tsol, Tres},
                           diffusion::AbstractDiffusion,
                           entropy_vars::AbstractVariables,
                           shockmesh::ShockedElements,
                           assembler::AssembleElementData) where {Tsol, Tres}

  # computeGradW computes Lambda * D * w, so all that remains to do is
  # compute Qx * grad_q_x
  # Note that this term is not entropy stable by itself, because Qx was
  # not replaced by -Qx^T + Ex.  The entire discretization should be
  # entropy-stable however.
  work = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.dim)
  Dx = zeros(mesh.numNodesPerElement, mesh.numNodesPerElement, mesh.dim)
  w_dot = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerElement)
  t1_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode, mesh.dim,
                 mesh.numNodesPerElement, mesh.numNodesPerElement)
  t2_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode, mesh.dim,
                 mesh.numNodesPerElement, mesh.numNodesPerElement)
  gradq_i = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement, mesh.dim)
  # if dim was the last dimension of t1_dot, we could re-use it, but instead
  # we have to use a 3rd array
  res_jac = eqn.params.calc_volume_integrals_data.res_jac
  op = SummationByParts.Subtract()

  for i=1:shockmesh.numShock
    i_full = shockmesh.elnums_all[i]

    q_i = ro_sview(eqn.q, :, :, i_full)
    w_i = ro_sview(capture.w_el, :, :, i)
    coords_i = ro_sview(mesh.coords, :, :, i_full)
    dxidx_i = ro_sview(mesh.dxidx, :, :, :, i_full)
    jac_i = ro_sview(mesh.jac, :, i_full)

    # to compute the jacobian, dq/dq = I, so dw/dq = the dw/dq evaluated at
    # each node
    for j=1:mesh.numNodesPerElement
      q_j = ro_sview(eqn.q, :, j, i_full)
      w_dot_j = sview(w_dot, :, :, j)
      getA0inv(entropy_vars, eqn.params, q_j, w_dot_j)
    end

    # apply Dx * w_dot
    calcDx(sbp, dxidx_i, jac_i, Dx)
    applyOperatorJac(Dx, w_dot, t1_dot)

    # need Dx * w for applyDiffusionTensor
    for d=1:mesh.dim
      Dx_d = ro_sview(Dx, :, :, d)
      gradq_d = sview(gradq_i, :, :, d)
      smallmatmatT!(w_i, Dx_d, gradq_d)
    end

    # apply the diffusion tensor to all nodes of the element
    applyDiffusionTensor_diff(diffusion, sbp, eqn.params, q_i, w_i, coords_i,
                              dxidx_i, jac_i, i, gradq_i, t1_dot, t2_dot)
#=
    # apply Qx and sum
    if eqn.params.use_Minv != 1
      @simd for d=1:mesh.dim
        @simd for j=1:mesh.numNodesPerElement
          @simd for i=1:mesh.numNodesPerElement
            Dx[i, j, d] *= sbp.w[i]/jac_i[i]  # Dx is now Qx
          end
        end
      end
    end
    applyOperatorJac(Dx, t2_dot, res_jac)
=#
    calcQxTransposed(sbp, dxidx_i, Dx)
    applyOperatorJac(Dx, t2_dot, res_jac, true, op)

    assembleElement(assembler, mesh, i_full, res_jac)
  end

  return nothing
end

function computeFaceTerm_diff(mesh, sbp, eqn, opts,
                      capture::SBPParabolicReduced2SC{Tsol, Tres},
                      shockmesh::ShockedElements, diffusion::AbstractDiffusion,
                      entropy_vars::AbstractVariables,
                      penalty::AbstractDiffusionPenalty,
                      assem::AssembleElementData) where {Tsol, Tres}

  delta_w = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  delta_w_dotL = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode,
                             mesh.numNodesPerFace, mesh.numNodesPerElement)
  delta_w_dotR = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode,
                             mesh.numNodesPerFace, mesh.numNodesPerElement)

  theta = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  theta_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  theta_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)

  res1L_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  res1L_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)

  res1R_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  res1R_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
#=
  res2L_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  res2L_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)

  res2R_dotL = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
  res2R_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                           mesh.numNodesPerFace, mesh.numNodesPerElement)
=#

  data = eqn.params.calc_face_integrals_data
  @unpack data res_jacLL res_jacLR res_jacRL res_jacRR

  #t1L = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  #t1R = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  t2L = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  t2R = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  subtract = SummationByParts.Subtract()
  add = SummationByParts.Add()

  for i=1:shockmesh.numInterfaces
    iface_red = shockmesh.ifaces[i].iface
    iface_idx = shockmesh.ifaces[i].idx_orig
    elnumL = shockmesh.elnums_all[iface_red.elementL]
    elnumR = shockmesh.elnums_all[iface_red.elementR]

    # get data needed for next steps
    wL = ro_sview(capture.w_el, :, :, iface_red.elementL)
    wR = ro_sview(capture.w_el, :, :, iface_red.elementR)
    qL = ro_sview(eqn.q, :, :, elnumL)
    qR = ro_sview(eqn.q, :, :, elnumR)
    gradwL = ro_sview(capture.grad_w, :, :, :, iface_red.elementL)
    gradwR = ro_sview(capture.grad_w, :, :, :, iface_red.elementR)

    coordsL = ro_sview(mesh.coords, :, :, elnumL)
    coordsR = ro_sview(mesh.coords, :, :, elnumR)
    nrm_face = ro_sview(mesh.nrm_face, :, :, iface_idx)
    dxidxL = ro_sview(mesh.dxidx, :, :, :, elnumL)
    dxidxR = ro_sview(mesh.dxidx, :, :, :, elnumR)
    jacL = ro_sview(mesh.jac, :, elnumL)
    jacR = ro_sview(mesh.jac, :, elnumR)
    resL = sview(eqn.res, :, :, elnumL)
    resR = sview(eqn.res, :, :, elnumR)

    # compute delta w tilde and theta_bar = Dgk w_k + Dgn w_n
    getFaceVariables_diff(eqn.params, capture, diffusion, entropy_vars,
                          sbp, mesh.sbpface,
                          iface_red, wL, wR, qL, qR, gradwL, gradwR,
                          coordsL, coordsR,
                          nrm_face, dxidxL, dxidxR, jacL, jacR,
                          delta_w, delta_w_dotL, delta_w_dotR)

    # apply the penalty coefficient matrix
    has_T4 = applyReduced2Penalty_diff(penalty, sbp, eqn.params,  mesh.sbpface,
                      diffusion,
                      iface_red, delta_w, delta_w_dotL, delta_w_dotR,
                      theta, theta_dotL, theta_dotR,
                      qL, qR, wL, wR, coordsL, coordsR,  nrm_face,
                      dxidxL, dxidxR, jacL, jacR,
                      res1L_dotL, res1L_dotR, res1R_dotL, res1R_dotR)

    # apply Rgk^T, Rgn^T, Dgk^T, Dgn^T
    fill!(res_jacLL, 0); fill!(res_jacRL, 0)
    fill!(res_jacLR, 0); fill!(res_jacRR, 0)
    interiorFaceIntegrate_jac!(mesh.sbpface, iface_red, res1L_dotL, res1R_dotL,
                               res_jacLL, res_jacRL, subtract, subtract, false)
    interiorFaceIntegrate_jac!(mesh.sbpface, iface_red, res1L_dotR, res1R_dotR,
                               res_jacLR, res_jacRR, subtract, subtract, false)

    iface_full = replace_interface(iface_red, elnumL, elnumR)
    if !has_T4
      assembleInterfaceVisc(assem, mesh.sbpface, mesh, iface_full, res_jacLL,
                            res_jacLR, res_jacRL, res_jacRR)
    else
      #TODO: the sparsity pattern of the Jacobian may be incorrect for this case
      assembleInterfaceFull(assem, mesh, iface_full, res_jacLL,
                            res_jacLR, res_jacRL, res_jacRR)
    end

  end  # end loop i

  return nothing
end

"""
  Differentiated version of [`getFaceVariables`](@ref).
"""
function getFaceVariables_diff(params::ParamType,
                          capture::SBPParabolicReduced2SC{Tsol, Tres},
                          diffusion::AbstractDiffusion,
                          entropy_vars::AbstractVariables,
                          sbp::AbstractOperator,
                          sbpface::AbstractFace, iface_red::Interface,
                          wL::AbstractMatrix, wR::AbstractMatrix,
                          qL::AbstractMatrix, qR::AbstractMatrix,
                          gradwL::Abstract3DArray, gradwR::Abstract3DArray,
                          coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                          nrm_face::AbstractMatrix,
                          dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                          jacL::AbstractVector, jacR::AbstractVector,
                          delta_w::AbstractMatrix,
                          delta_w_dotL::Abstract4DArray, delta_w_dotR::Abstract4DArray,
                         ) where {Tsol, Tres}

  numDofPerNode, numNodesPerElement = size(wL)
  numNodesPerFace = size(delta_w, 2)
  dim = size(gradwL, 3)

  getFaceVariables(capture, sbpface, iface_red, wL, wR, gradwL, gradwR,
                   nrm_face, delta_w)

  @unpack capture wL_dot wR_dot Dx t1 t1_dot t2L_dot t2R_dot t3L_dot t3R_dot

  # derivative of delta_w wrt qL and qR
  for i=1:numNodesPerElement
    wL_dot_i = sview(wL_dot, :, :, i)
    wR_dot_i = sview(wR_dot, :, :, i)
    qL_i = ro_sview(qL, :, i)
    qR_i = ro_sview(qR, :, i)
    getA0inv(entropy_vars, params, qL_i, wL_dot_i)
    getA0inv(entropy_vars, params, qR_i, wR_dot_i)
  end

  fill!(delta_w_dotL, 0); fill!(delta_w_dotR, 0)
  interiorFaceInterpolate_jac!(sbpface, iface_red, wL_dot, wR_dot,
                              delta_w_dotL, delta_w_dotR)

  @simd for i=1:length(delta_w_dotR)
    delta_w_dotR[i] = -delta_w_dotR[i]  # wR is negated in wL - wR
  end

  return nothing
end

#------------------------------------------------------------------------------
# Dirichlet boundary condition

"""
  Differentiated Dirichlet boundary term
"""
function computeDirichletBoundaryTerm_diff(mesh, sbp, eqn, opts,
                      capture::ParabolicReduceds{Tsol, Tres},
                      shockmesh::ShockedElements,
                      assem::AssembleElementData
                      ) where {Tsol, Tres}

  for i=1:mesh.numBC
    bc_range = (shockmesh.bndry_offsets[i]:(shockmesh.bndry_offsets[i+1]-1))
    bc_func = mesh.bndry_funcs[i]

    calcBoundaryFlux_diff(mesh, sbp, eqn, opts, shockmesh, capture,
                     capture.penalty, capture.diffusion,
                     capture.entropy_vars, bc_func, bc_range, assem)
  end

  return nothing
end

function calcBoundaryFlux_diff(mesh::AbstractMesh, sbp, eqn::EulerData, opts,
                      shockmesh::ShockedElements,
                      capture::ParabolicReduceds{Tsol, Tres},
                      penalty::AbstractDiffusionPenalty,
                      diffusion::AbstractDiffusion,
                      entropy_vars::AbstractVariables, bc_func::BCType,
                      bc_range::AbstractVector,
                      assem::AssembleElementData,
                      ) where {Tsol, Tres}

  delta_w = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  delta_w_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                            mesh.numNodesPerFace, mesh.numNodesPerElement)
  delta_w_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                             mesh.numNodesPerFace, mesh.numNodesPerElement)
  res1 = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  res1_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                         mesh.numNodesPerFace, mesh.numNodesPerElement)

  res_dot = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                        mesh.numNodesPerElement, mesh.numNodesPerElement)
  res_dotR = zeros(Tres, mesh.numDofPerNode, mesh.numDofPerNode,
                        mesh.numNodesPerElement, mesh.numNodesPerElement)

  #op = SummationByParts.Subtract()
  subtract = SummationByParts.Subtract()
  op = SummationByParts.Add()

  for i in bc_range
    bndry_i = shockmesh.bndryfaces[i].bndry
    idx_orig = shockmesh.bndryfaces[i].idx_orig
    elnum_orig = shockmesh.elnums_all[bndry_i.element]
    fill!(res_dot, 0)

    w_i = ro_sview(capture.w_el, :, :, bndry_i.element)
    q_i = ro_sview(eqn.q, :, :, elnum_orig)
    coords_i = ro_sview(mesh.coords_bndry, :, :, idx_orig)
    nrm_i = ro_sview(mesh.nrm_bndry, :, :, idx_orig)
    #alpha = capture.alpha_b[i]
    alpha = 1
    dxidxL = ro_sview(mesh.dxidx, :, :, :, elnum_orig)
    jacL = ro_sview(mesh.jac, :, elnum_orig)
    res_i = sview(eqn.res, :, :, elnum_orig)

    computeDirichletDelta_diff(capture, eqn.params, eqn.params_complex,
                          mesh.sbpface, bndry_i,
                          bc_func, entropy_vars, w_i, q_i, coords_i,
                          nrm_i, delta_w, delta_w_dot)

    applyDirichletPenalty_diff(penalty, sbp, eqn.params, mesh.sbpface, diffusion,
                          bndry_i, delta_w, delta_w_dot, q_i, w_i, coords_i,
                          nrm_i, alpha,
                          dxidxL, jacL, res1_dot)

    # apply R^T to T_D * delta_w
    #fill!(res_dot, 0)
    boundaryFaceIntegrate_jac!(mesh.sbpface, bndry_i.face, res1_dot, res_dot,
                               subtract, include_quadrature=false)
#    boundaryFaceIntegrate_jac!(mesh.sbpface, bndry_i.face, delta_w_dot, res_dot,
#                               op, include_quadrature=false)

    # apply B and then Dgk^T to delta_w
    for q=1:mesh.numNodesPerElement
      for p=1:mesh.numNodesPerFace
        for j=1:mesh.numDofPerNode
          for i=1:mesh.numDofPerNode
            delta_w_dot[i, j, p, q] *= -mesh.sbpface.wface[p]
          end
        end
      end
    end

    for p=1:mesh.numNodesPerFace
      for i=1:mesh.numDofPerNode
        delta_w[i,p] *= -mesh.sbpface.wface[p]
      end
    end

    # the dotR variables aren't needed for boundary terms
    applyDgkTranspose_diff(capture, sbp, eqn.params, mesh.sbpface, bndry_i,
                      diffusion,
                      delta_w, delta_w_dot, delta_w_dotR, q_i, w_i, coords_i,
                      nrm_i, dxidxL, jacL, res_dot, res_dotR, subtract)

    assembleElement(assem, mesh, elnum_orig, res_dot)
  end

  return nothing
end

function computeDirichletDelta_diff(capture::ParabolicReduceds{Tsol, Tres},
                               params::ParamType, params_c::ParamType,
                               sbpface::AbstractFace, bndry::Boundary,
                               func::BCType, entropy_vars::AbstractVariables,
                               wL::AbstractMatrix, qL::AbstractMatrix,
                               coords::AbstractMatrix, nrm_face::AbstractMatrix,
                               delta_w::AbstractMatrix,
                               delta_w_dot::Abstract4DArray
                              ) where {Tsol, Tres}

  numDofPerNode, numNodesPerFace = size(delta_w)
  numNodesPerElement = size(wL, 2)

  #TODO: maybe interpolating the conservative variables and then converting
  #      would be more consistent, but Eq. 31 requires the same interpolation
  #      Rgk * w for both the interface and the boundary term.
  #      For the boundary state qg we have more leeway because it is zero
  #      for the energy stability analysis (-> entropy stability when u = w)
  wface = zeros(Tsol, numDofPerNode, numNodesPerFace)
  boundaryFaceInterpolate!(sbpface, bndry.face, wL, wface)

  # jacobian
  wL_dot = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerElement)
#  wface_dot = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerFace,
#                    numNodesPerElement)
  for i=1:numNodesPerElement
    wL_dot_i = sview(wL_dot, :, :, i)
    qL_i = ro_sview(qL, :, i)
    getA0inv(entropy_vars, params, qL_i, wL_dot_i)
  end
#  println("qL = \n", real(qL))
#  println("wL_dot = \n", real(wL_dot))
#  println("perm = ", sbpface.perm[:, bndry.face])
  fill!(delta_w_dot, 0)
  boundaryFaceInterpolate_jac!(sbpface, bndry.face, wL_dot, delta_w_dot)

  qface = zeros(Tsol, numDofPerNode, numNodesPerFace)
  qgface = zeros(Tres, numDofPerNode, numNodesPerFace)
  wgface = zeros(Tres, numDofPerNode)
  aux_vars = Tres[]

  # the inviscid boundary conditions interpolate q to the face and then
  # compute qg, so do that here as well for consistency
  boundaryFaceInterpolate!(sbpface, bndry.face, qL, qface)
  for i=1:numNodesPerFace
    q_i = sview(qface, :, i)
    qgface_i = sview(qgface, :, i)  # need this later
    coords_i = sview(coords, :, i)
    nrm_i = sview(nrm_face, :, i)

    getDirichletState(func, params, q_i, aux_vars, coords_i, nrm_i, qgface_i)

    convertToEntropy(entropy_vars, params, qgface_i, wgface)

    # compute the delta
    for j=1:numDofPerNode
      delta_w[j, i] = wface[j, i] - wgface[j]
    end
  end

  # jacobian
  # Compute dqg/qface first, then R^T * dqg/qface to get dqg/dqvolume

  op = SummationByParts.Subtract()
  qg_dot = zeros(Tsol, numDofPerNode, numDofPerNode)
  A0inv = zeros(Tsol, numDofPerNode, numDofPerNode)
  wg_dotface = zeros(Tsol, numDofPerNode, numDofPerNode, numNodesPerFace)
  q_c = zeros(Complex128, numDofPerNode)
  qg_c = zeros(Complex128, numDofPerNode)
  h = 1e-20
  pert = Complex128(0, h)
  for i=1:numNodesPerFace
    for j=1:numDofPerNode
      q_c[j] = qface[j, i]
    end
    coords_i = sview(coords, :, i)
    nrm_i = sview(nrm_face, :, i)

    # compute jacobian of qgface wrt qface
    for j=1:numDofPerNode
      q_c[j] += pert

      getDirichletState(func, params_c, q_c, aux_vars, coords_i, nrm_i, qg_c)
      for k=1:numDofPerNode
        qg_dot[k, j] = imag(qg_c[k])/h
      end

      q_c[j] -= pert
    end

    # jacobian of convertToEntropy
    getA0inv(entropy_vars, params, sview(qgface, :, i), A0inv)
    wg_dot_i = sview(wg_dotface, :, :, i)
    smallmatmat!(A0inv, qg_dot, wg_dot_i)
  end

  # apply R^T to wg_dotface to get d wg/dqL, and subtract that from delta_w_dot
  # interpolation is linear, so apply the regular operation repatedly
  # It would be better if SBP had a jacobian operation for this
  faceToVolume_jac!(sbpface, bndry.face, wg_dotface, delta_w_dot, op)

  return nothing
end

"""
  Jacobian of Dgk * t2L only
"""
function applyDgkTranspose_diff(capture::ParabolicReduceds{Tsol, Tres}, sbp,
                       params::ParamType,
                       sbpface, iface::Union{Interface, Boundary},
                       diffusion::AbstractDiffusion,
                       t2L::AbstractMatrix, t2L_dotL::Abstract4DArray,
                       t2L_dotR::Abstract4DArray,
                       qL::AbstractMatrix,
                       wL::AbstractMatrix,
                       coordsL::AbstractMatrix,
                       nrm_face::AbstractMatrix,
                       dxidxL::Abstract3DArray,
                       jacL::AbstractVector,
                       resL_dotL::Abstract4DArray, resL_dotR::Abstract4DArray,
                       op::SummationByParts.UnaryFunctor=SummationByParts.Add()
                      ) where {Tsol, Tres}

  dim, numNodesPerFace = size(nrm_face)
  numNodesPerElement = size(dxidxL, 3)
  numDofPerNode = size(wL, 1)
  add = SummationByParts.Add()
  subtract = SummationByParts.Subtract()

  @unpack capture temp1L temp2L
  fill!(temp2L, 0);
  @unpack capture t3L_dotL t3L_dotR Dx
  @unpack capture t4L_dotL t4L_dotR
  @unpack capture t5L_dotL t5L_dotR

  # apply N and R^T
  @simd for d=1:dim
    @simd for j=1:numNodesPerFace
      @simd for k=1:numDofPerNode
        temp1L[k, j] =  nrm_face[d, j]*t2L[k, j]
      end
    end

    tmp2L = sview(temp2L, :, :, d)
    boundaryFaceInterpolate_rev!(sbpface, getFaceL(iface), tmp2L, temp1L)
  end

  # multiply by N
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerFace
      @simd for d=1:dim
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            t3L_dotL[i, j, d, p, q] = nrm_face[d, p]*t2L_dotL[i, j, p, q]
            t3L_dotR[i, j, d, p, q] = nrm_face[d, p]*t2L_dotR[i, j, p, q]
          end
        end
      end
    end
  end

  # apply R^T
  fill!(t4L_dotL, 0); fill!(t4L_dotR, 0)
  #TODO: get rid of keyword arguments
  boundaryFaceIntegrate_jac!(sbpface, getFaceL(iface), t3L_dotL, t4L_dotL,
                             include_quadrature=false)
  boundaryFaceIntegrate_jac!(sbpface, getFaceL(iface), t3L_dotR, t4L_dotR,
                             include_quadrature=false)

  # multiply by D^T Lambda
  applyDiffusionTensor_diff(diffusion, sbp, params, qL, wL, coordsL, dxidxL,
                            jacL, getElementL(iface), temp2L,
                            t4L_dotL, t4L_dotR, t5L_dotL, t5L_dotR)

  calcDxTransposed(sbp, dxidxL, jacL, Dx)
  applyOperatorJac(Dx, t5L_dotL, resL_dotL, false, op)
  applyOperatorJac(Dx, t5L_dotR, resL_dotR, false, op)

  return nothing
end

#------------------------------------------------------------------------------
# Penalty

"""
  Differentiated version of BR2 penalty
"""
function applyReduced2Penalty_diff(penalty::BR2Penalty{Tsol, Tres}, sbp,
                      params::AbstractParamType, sbpface,
                      diffusion::AbstractDiffusion, iface::Interface,
                      delta_w::AbstractMatrix{Tsol},
                      delta_w_dotL::Abstract4DArray, delta_w_dotR::Abstract4DArray,
                      theta::AbstractMatrix{Tres},
                      theta_dotL::Abstract4DArray, theta_dotR::Abstract4DArray,
                      qL_el::AbstractMatrix, qR_el::AbstractMatrix,
                      wL::AbstractMatrix, wR::AbstractMatrix,
                      coordsL::AbstractMatrix, coordsR::AbstractMatrix,
                      nrm_face::AbstractMatrix,
                      dxidxL::Abstract3DArray, dxidxR::Abstract3DArray,
                      jacL::AbstractVector, jacR::AbstractVector,
                      res1L_dotL::Abstract4DArray, res1L_dotR::Abstract4DArray,
                      res1R_dotL::Abstract4DArray, res1R_dotR::Abstract4DArray,
                     ) where {Tsol, Tres}

  numDofPerNode, numNodesPerFace = size(delta_w)
  numNodesPerElement = length(jacL)
  dim = size(nrm_face, 1)
  add = SummationByParts.Add()

  # nodes in the stencil of R
  nodesL = sview(sbpface.perm, :, iface.faceL)
  nodesR = sview(sbpface.perm, :, iface.faceR)

  # the first L/R indicates if this is element kappa or nu
  # the second L/R indices if the derivative is wrt qL or qR
  @unpack penalty delta_w_n qL qR t1_dotL t1_dotR t2LL t2LR t2RL t2RR
  @unpack penalty t3LL t3LR t3RL t3RR t4RL t4RR
  t4LL = t1_dotL; t4LR = t1_dotR

  fill!(qL, 0); fill!(qR, 0)

  #--------------------------
  # T1

  # compute primal quantity needed later
  # multiply by normal vector, then R^T B
  #alpha_g = 1/(dim + 1)  # = 1/number of faces of a simplex
  @simd for d1=1:dim
    @simd for j=1:numNodesPerFace
      @simd for k=1:numDofPerNode
        delta_w_n[k, j] = 0.25*delta_w[k, j]*nrm_face[d1, j]
      end
    end

    qL_d = sview(qL, :, :, d1); qR_d = sview(qR, :, :, d1)
    interiorFaceIntegrate!(sbpface, iface, delta_w_n, qL_d, qR_d)

    # interiorFaceIntegrates -= the second output
    @simd for j=1:numNodesPerElement
      @simd for k=1:numDofPerNode
        qR_d[k, j] = -qR_d[k, j]
      end
    end
  end

  # compute derivative quantity
  # apply B and [Nx, Ny]
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerFace
      @simd for d=1:dim
        fac = 0.25*sbpface.wface[p]*nrm_face[d, p]
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            t1_dotL[i, j, d, p, q] = fac*delta_w_dotL[i, j, p, q]
            t1_dotR[i, j, d, p, q] = fac*delta_w_dotR[i, j, p, q]
          end
        end
      end
    end
  end

  # interpolate back to volume nodes
  fill!(t2LL, 0); fill!(t2RL, 0); fill!(t2LR, 0); fill!(t2RR, 0)
  interiorFaceIntegrate_jac!(sbpface, iface, t1_dotL, t1_dotL, t2LL, t2RL, add,
                             add, false)
  interiorFaceIntegrate_jac!(sbpface, iface, t1_dotR, t1_dotR, t2LR, t2RR, add,
                             add, false)

  # apply diffusion tensor
  applyDiffusionTensor_diff(diffusion, sbp, params, qL_el, wL, coordsL, dxidxL,
                            jacL, iface.elementL, sbpface,
                            iface.faceL, qL, t2LL, t2LR, t3LL, t3LR)
  # reverse RR and RL here because Lambda is a function of qR, not qL
  applyDiffusionTensor_diff(diffusion, sbp, params, qR_el, wR, coordsR, dxidxR,
                            jacR, iface.elementR, sbpface,
                            iface.faceR, qR, t2RR, t2RL, t3RR, t3RL)

  # apply inverse mass matrix (only nodes needed for interpolation later)
  @simd for q=1:numNodesPerElement
    # elementL
    @simd for p in nodesL
      facL = jacL[p]/sbp.w[p]
      @simd for d=1:dim
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            t3LL[i, j, d, p, q] *= facL
            t3LR[i, j, d, p, q] *= facL
          end
        end
      end
    end
    # elementR
    @simd for p in nodesR
      facR = jacR[p]/sbp.w[p]
      @simd for d=1:dim
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            t3RL[i, j, d, p, q] *= facR
            t3RR[i, j, d, p, q] *= facR
          end
        end
      end
    end
  end  # end q

  # interpolate back to the face
  fill!(t4LL, 0); fill!(t4RL, 0); fill!(t4LR, 0); fill!(t4RR, 0)
  interiorFaceInterpolate_jac!(sbpface, iface, t3LL, t3RL, t4LL, t4RL)
  interiorFaceInterpolate_jac!(sbpface, iface, t3LR, t3RR, t4LR, t4RR)

  # apply B * [Nx, Ny] and sum
#  fill!(res1L_dotL, 0); fill!(res1L_dotR, 0)  #TODO: tile this
#  fill!(res1R_dotL, 0); fill!(res1R_dotR, 0)
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerFace
      @simd for j=1:numDofPerNode
        @simd for i=1:numDofPerNode
          res1L_dotL[i, j, p, q] = 0
          res1L_dotR[i, j, p, q] = 0
          res1R_dotL[i, j, p, q] = 0
          res1R_dotR[i, j, p, q] = 0
        end
      end

      @simd for d=1:dim
        fac = sbpface.wface[p]*nrm_face[d, p]
        @simd for j=1:numDofPerNode
          @simd for i=1:numDofPerNode
            valL = fac*(t4LL[i, j, d, p, q] + t4RL[i, j, d, p, q])
            valR = fac*(t4LR[i, j, d, p, q] + t4RR[i, j, d, p, q])
            #println("t2L_d_dot = ", real(t4LL[i, j, d, p, q]))
            #println("t2R_d_dot = ", real(t4RL[i, j, d, p, q]))
            #println("val = ", real(valL))
            res1L_dotL[i, j, p, q] +=  valL
            res1L_dotR[i, j, p, q] +=  valR
            res1R_dotL[i, j, p, q] += -valL  # delta_w is reversed for elementR
            res1R_dotR[i, j, p, q] += -valR
          end
        end
      end
    end
  end

  return false
end


text:
using Random
Random.seed!(209)

@testset "0D grids" begin
  @test Space((1,)).graph == Agents.Graph(1)
end

@testset "1D grids" begin
  a = Space((5,1))
  ae = collect(Agents.LightGraphs.edges(a.graph))

  b = Space((5,1), periodic=true)
  be = collect(Agents.LightGraphs.edges(b.graph))

  @test ae == [Agents.LightGraphs.Edge(1,2),Agents.LightGraphs.Edge(2,3), Agents.LightGraphs.Edge(3,4), Agents.LightGraphs.Edge(4,5)]
  @test be == [Agents.LightGraphs.Edge(1,2), Agents.LightGraphs.Edge(1,5), Agents.LightGraphs.Edge(2,3), Agents.LightGraphs.Edge(3,4), Agents.LightGraphs.Edge(4,5)]
end

@testset "2D grids" begin
  @test Space((2,1)).graph == Agents.Graph(2,1)

  a = Space((2,3))
  b = Space((2,3), periodic=true) # 2D grid

  @test a.dimensions == (2,3)
  @test b.dimensions == (2,3)
  @test length(a.agent_positions) == 6
  @test length(b.agent_positions) == 6

  @test Agents.nv(a) == 6
  @test Agents.ne(a) == 7
  @test Agents.nv(b) == 6
  @test Agents.ne(b) == 9

  ae = collect(Agents.LightGraphs.edges(a.graph))
  be = collect(Agents.LightGraphs.edges(b.graph))

  @test ae == [Agents.LightGraphs.Edge(1,2), Agents.LightGraphs.Edge(1,3), Agents.LightGraphs.Edge(2,4), Agents.LightGraphs.Edge(3,4), Agents.LightGraphs.Edge(3,5), Agents.LightGraphs.Edge(4,6), Agents.LightGraphs.Edge(5,6)]
  @test be == [Agents.LightGraphs.Edge(1,2), Agents.LightGraphs.Edge(1,3), Agents.LightGraphs.Edge(1,5), Agents.LightGraphs.Edge(2,4), Agents.LightGraphs.Edge(2,6), Agents.LightGraphs.Edge(3,4), Agents.LightGraphs.Edge(3,5), Agents.LightGraphs.Edge(4,6), Agents.LightGraphs.Edge(5,6)]
end

@testset "2D Moore" begin
  a = Space((3, 3), moore=true)
  @test Agents.nv(a) == 9
  @test Agents.ne(a) == 20

  a = Space((3, 4), moore=true)
  @test Agents.nv(a) == 12
  @test Agents.ne(a) == 29

  b = Space((3, 2), moore=true, periodic=true)
  @test Agents.nv(b) == 6
  @test Agents.ne(b) == 15

  b = Space((3, 3), moore=true, periodic=true)
  @test Agents.nv(b) == 9
  @test Agents.ne(b) == 36
end

@testset "3D grid" begin
  g1 = Space((2,3,2))
  g2 = Space((2,3,3))
  g3 = Space((2,3,2), periodic=true)
  g4 = Space((2,3,3), periodic=true)
  @test Agents.ne(g1) == 20
  @test Agents.ne(g2) == 33
  @test Agents.ne(g3) == 24
  @test Agents.ne(g4) == 39
end

@testset "grid coord/vertex conversions" begin
  @test coord2vertex((2,2,1), (3,4,1)) == 5
  @test coord2vertex((2,2), (3,4)) == 5
  @test coord2vertex((2,2,1), (3,4,1)) == 5
  @test coord2vertex((2,2), (3,4)) == 5
  @test vertex2coord(5, (3,4,1)) == (2,2,1)
  @test vertex2coord(5, (3,4)) == (2,2)

  @test coord2vertex((2,2,1), (2,3,3)) == 4
  @test coord2vertex((2,2,2), (2,3,3)) == 10
  @test coord2vertex((2,2,3), (2,3,3)) == 16
  @test coord2vertex(((1,3,1)), (2,3,3)) == 5
  @test coord2vertex((1,3,2), (2,3,3)) == 11
  @test coord2vertex((1,3,3), (2,3,3)) == 17

  @test vertex2coord(5, (2,3,3)) == (1,3,1)
  @test vertex2coord(7, (2,3,3)) == (1,1,2)
  @test vertex2coord(13, (2,3,3)) == (1,1,3)
  @test vertex2coord(13, (2,3,3)) == (1,1,3)
  @test vertex2coord(15, (2,3,3)) == (1,2,3)
  @test vertex2coord(18, (2,3,3)) == (2,3,3)
  @test vertex2coord(18, (2,3,3)) == (2,3,3)
end

@testset "nodes" begin
  space = Space((3,3))
  model = ABM(Agent1, space)
  for node in nodes(model)
    if rand() > 0.7
      add_agent!(node, model)
    end
  end
  @test nodes(model) == 1:9
  @test nodes(model, by=:random) == [8, 9, 1, 7, 6, 5, 3, 2, 4]
  @test nodes(model, by=:population) == [4, 5, 8, 9, 1, 2, 3, 6, 7]
  @test length(get_node_contents(4, model)) > length(get_node_contents(7, model))
  @test_throws ErrorException nodes(model, by=:notreal)
end


text:
# Tests for the total correlation function

total_correlation = entropy1 + entropy2 + entropy3 - entropy123

# Test total correlation
@test get_total_correlation(arr1, arr2, arr3) ≈ total_correlation
println("Total correlation passed.")

# Test total correlation formula
@test apply_total_correlation_formula(entropy1, entropy2, entropy3, entropy123) ≈ total_correlation
println("Total correlation formula passed.")


text:
using BinaryBuilder, Pkg
using BinaryBuilder: download_source

# Read in input `.json` file
json = String(read(ARGS[1]))
buff = IOBuffer(strip(json))
objs = []
while !eof(buff)
    push!(objs, BinaryBuilder.JSON.parse(buff))
end

# Merge the multiple outputs into one
merged = BinaryBuilder.merge_json_objects(objs)
BinaryBuilder.cleanup_merged_object!(merged)

# Download all sources, unless we're in a skip build situation
if get(ENV, "SKIP_BUILD", "false") != "true"
    download_source.(merged["sources"]; verbose=true)
end

# Also initialize JLL package directories
src_name = merged["name"]
code_dir = joinpath(Pkg.devdir(), "$(src_name)_jll")
# Always start from a clean slate
rm(code_dir; recursive=true, force=true)

# Then export platforms to file
open(ARGS[2], "w") do io
    println(io, join(triplet.(merged["platforms"]), " "))
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
export AbstractParameter, Parameter, value

abstract type AbstractParameter{T} end
const AP{T} = AbstractParameter{T}

# Sometime when setup our models, we require it's parameters to follow
# certain constrains, e.g. the support for scaling factor `σ` & length scale `l`
# of a GP kernel should be ℜ⁺, however, adding those constrains make the optimization
# process less convenient. `Parameter` is used to handle this, it's bijector is used
# to map constrained variable to unconstrained one, which is used during the optimization,
# and the inverse of the bijector does the reverse, which will provide us with the constrained
# variables we need.
struct Parameter{T, N, fT<:Bijector{N}} <: AP{T}
    x::AbstractArray{T}
    f::fT

    function Parameter{T}(y, f::fT) where {T<:Real, fT}
        N = ndims(y)
        x = inv(f)(y)
        if N == 0
            return new{T, N, fT}(T[x], f)
        else
            return new{T, N, fT}(x, f)
        end
    end
end

Parameter(y, f) = Parameter{eltype(y)}(y, f)

# interface
value(y::Parameter{T, 0}) where {T} = y.f(first(y.x))
value(y::Parameter) = y.f(y.x)

# if no constrain is added, it will use `Identity` bijector by default
Parameter(y) = Parameter(y, Identity{ndims(y)}())

# if we want positivity constrain, then use `Exp` bijector
Parameter(y, ::Val{:pos}) = Parameter(y, Bijectors.Exp{ndims(y)}())


text:
@testset "376.wiggle-subsequence.jl" begin
    @test wiggle_max_length([1,2,3,4,5,6,7,8,9]) == 2
    @test wiggle_max_length([1,17,5,10,13,15,10,5,16,8]) == 7
    @test wiggle_max_length([1,7,4,9,2,5]) == 6
    @test wiggle_max_length([0, 0]) == 1
end


text:
using tempfile: NamedTemporaryFile
if abspath(PROGRAM_FILE) == @__FILE__
    NamedTempFile::new() do temp_file
        file_path = temp_file.name
        open(file_path, "w") do f
            write(f, "hello")
        end
        open(file_path, "r") do f
            @assert(read(f, 1) == "h")
            @assert(read(f) == "ello")
            println("OK")
        end
    end
end


text:
using FrameConnection
using Test

@testset "FrameConnection.jl" begin
    # Write your tests here.
end


text:
if VERSION < v"0.5.0"
    # JuliaLang/julia#16219
    if !isdefined(Base, Symbol("@static"))
        macro static(ex)
            if isa(ex, Expr)
                if ex.head === :if
                    cond = eval(current_module(), ex.args[1])
                    if cond
                        return esc(ex.args[2])
                    elseif length(ex.args) == 3
                        return esc(ex.args[3])
                    else
                        return nothing
                    end
                end
            end
            throw(ArgumentError("invalid @static macro"))
        end
    end
    @eval is_windows() = $(OS_NAME == :Windows)
    @eval is_apple()   = $(OS_NAME == :Darwin)
    @eval is_linux()   = $(OS_NAME == :Linux)
    @eval is_bsd()     = $(OS_NAME in (:FreeBSD, :OpenBSD, :NetBSD, :Darwin, :Apple))
    @eval is_unix()    = $(is_linux() || is_bsd())

    const ASCIIStr  = ASCIIString
    const UTF8Str   = UTF8String
    const UniStr    = UTF16String
    const ByteStr   = ByteString

    const cvt_utf8  = utf8
    const cvt_utf16 = utf16
else
    using LegacyStrings
    const ASCIIStr  = LegacyStrings.ASCIIString
    const UTF8Str   = LegacyStrings.UTF8String
    const UniStr    = LegacyStrings.UTF16String
    const ByteStr   = Union{ASCIIStr, UTF8Str, String}

    const cvt_utf8  = LegacyStrings.utf8
    const cvt_utf16 = LegacyStrings.utf16
end


text:
module WordNet

using Compat

include("pointer.jl")
include("lemma.jl")
include("constants.jl")
include("synset.jl")
include("db.jl")
include("operations.jl")
include("sensekeys.jl")

end


text:
struct XVG
    Data::Array{Array{Float64,1},1}
    x_label::SubString{String}
    y_label::SubString{String}
    plot_labels
    title::SubString{String}
end

function read_xvg(path::AbstractString)
    a = readlines(path)
    plot_atributes = filter(x->startswith(x,"@"),a)
    if length(plot_atributes) >= 11
            if length(plot_atributes) > 11
                plot_label = Array{AbstractString,1}(undef,length(plot_atributes)-11 + 1)
                for i in 11:length(plot_atributes)
                    plot_label[i-10] = replace(replace(replace(strip(split(plot_atributes[i],"legend")[2]),"\\s"=>"_"),"\""=>""),"\\N"=>"")
                end
                String.(plot_label)
            else
                plot_label = split(plot_atributes[11],"legend")[2][3:end-1]
            end
    else
        plot_label = split(plot_atributes[1],"title")[2][3:end-1]
    end
    title = split(plot_atributes[1],"title")[2][3:end-1]
    x_label = split(plot_atributes[2],"label")[2][3:end-1]
    y_label = split(plot_atributes[3],"label")[2][3:end-1]
    dp = map(x->parse.(Float64,x),split.(filter(x->startswith(x," "),a)))

    return XVG(dp,x_label,y_label,plot_label,title)
end

using RecipesBase

@recipe function f(a::XVG)
    label --> false
    xguide --> a.x_label
    yguide --> a.y_label
    title --> a.title
    legend --> :outertopright

    if typeof(a.plot_labels) != SubString{String}
        labels --> reshape(a.plot_labels,(1,length(a.plot_labels)))
    else
        labels --> a.plot_labels
    end
    @series begin
        if a.Data[1] |> length == 2
            getindex.(a.Data,1),getindex.(a.Data,2)
        else
            getindex.(a.Data,1),[getindex.(a.Data,2), getindex.(a.Data,3),getindex.(a.Data,4),getindex.(a.Data,5)]
        end

    end
end


text:
function HDA(g, threshold=sqrt(nv(g)))
	attack_nodes = high_degree_adaptive(g)
	present = fill(true, nv(g))
	present[attack_nodes] .= false
	reinsert_nodes = reverse_greedy!(present, g, threshold)
	return [i for i in attack_nodes if !present[i]]
end

function CoreHD(g, threshold=sqrt(nv(g)))
	decycling_nodes = decore(g)
	present = fill(true, nv(g))
	present[decycling_nodes] .= false
	treebreak_nodes, _ = tree_break!(present, g, threshold)
	reinsert_nodes = reverse_greedy!(present, g, threshold)
	return [i for i in [decycling_nodes;treebreak_nodes] if !present[i]]
end

function WeakNei(g, threshold=sqrt(nv(g)))
	decycling_nodes = weak_neighbor_decore(g)
	present = fill(true, nv(g))
	present[decycling_nodes] .= false
	treebreak_nodes, _ = tree_break!(present, g, threshold)
	reinsert_nodes = reverse_greedy!(present, g, threshold)
	return [i for i in [decycling_nodes;treebreak_nodes] if !present[i]]
end


text:
@testset "509.fibonacci-number.jl" begin
    @test fib(5) == 5
end


text:
using ITensors, HDF5, Test

include("util.jl")

@testset "HDF5 Read and Write" begin
  i = Index(2, "i")
  j = Index(3, "j")
  k = Index(4, "k")

  @testset "TagSet" begin
    ts = TagSet("A,Site,n=2")
    fo = h5open("data.h5", "w")
    write(fo, "tags", ts)
    close(fo)

    fi = h5open("data.h5", "r")
    rts = read(fi, "tags", TagSet)
    close(fi)
    @test rts == ts
  end

  @testset "Index" begin
    i = Index(3, "Site,S=1")
    fo = h5open("data.h5", "w")
    write(fo, "index", i)
    close(fo)

    fi = h5open("data.h5", "r")
    ri = read(fi, "index", Index)
    close(fi)
    @test ri == i

    # primed Index
    i = Index(3, "Site,S=1")
    i = prime(i, 2)
    fo = h5open("data.h5", "w")
    write(fo, "index", i)
    close(fo)

    fi = h5open("data.h5", "r")
    ri = read(fi, "index", Index)
    close(fi)
    @test ri == i
  end

  @testset "IndexSet" begin
    is = IndexSet(i, j, k)

    fo = h5open("data.h5", "w")
    write(fo, "inds", is)
    close(fo)

    fi = h5open("data.h5", "r")
    ris = read(fi, "inds", IndexSet)
    close(fi)
    @test ris == is
  end

  @testset "Dense ITensor" begin

    # default constructed case
    T = ITensor()

    fo = h5open("data.h5", "w")
    write(fo, "defaultT", T)
    close(fo)

    fi = h5open("data.h5", "r")
    rT = read(fi, "defaultT", ITensor)
    close(fi)
    @test typeof(storage(T)) == typeof(storage(ITensor()))

    # real case
    T = randomITensor(i, j, k)

    fo = h5open("data.h5", "w")
    write(fo, "T", T)
    close(fo)

    fi = h5open("data.h5", "r")
    rT = read(fi, "T", ITensor)
    close(fi)
    @test norm(rT - T) / norm(T) < 1E-10

    # complex case
    T = randomITensor(ComplexF64, i, j, k)

    fo = h5open("data.h5", "w")
    write(fo, "complexT", T)
    close(fo)

    fi = h5open("data.h5", "r")
    rT = read(fi, "complexT", ITensor)
    close(fi)
    @test norm(rT - T) / norm(T) < 1E-10
  end

  @testset "QN ITensor" begin
    i = Index(QN("A", -1) => 3, QN("A", 0) => 4, QN("A", +1) => 3; tags="i")
    j = Index(QN("A", -2) => 2, QN("A", 0) => 3, QN("A", +2) => 2; tags="j")
    k = Index(QN("A", -1) => 1, QN("A", 0) => 1, QN("A", +1) => 1; tags="k")

    # real case
    T = randomITensor(QN("A", 1), i, j, k)

    fo = h5open("data.h5", "w")
    write(fo, "T", T)
    close(fo)

    fi = h5open("data.h5", "r")
    rT = read(fi, "T", ITensor)
    close(fi)
    @test rT ≈ T

    # complex case
    T = randomITensor(ComplexF64, i, j, k)

    fo = h5open("data.h5", "w")
    write(fo, "complexT", T)
    close(fo)

    fi = h5open("data.h5", "r")
    rT = read(fi, "complexT", ITensor)
    close(fi)
    @test rT ≈ T
  end

  @testset "MPO/MPS" begin
    N = 6
    sites = siteinds("S=1/2", N)

    # MPO
    mpo = makeRandomMPO(sites)

    fo = h5open("data.h5", "w")
    write(fo, "mpo", mpo)
    close(fo)

    fi = h5open("data.h5", "r")
    rmpo = read(fi, "mpo", MPO)
    close(fi)
    @test prod([norm(rmpo[i] - mpo[i]) / norm(mpo[i]) < 1E-10 for i in 1:N])

    # MPS
    mps = makeRandomMPS(sites)
    fo = h5open("data.h5", "w")
    write(fo, "mps", mps)
    close(fo)

    fi = h5open("data.h5", "r")
    rmps = read(fi, "mps", MPS)
    close(fi)
    @test prod([norm(rmps[i] - mps[i]) / norm(mps[i]) < 1E-10 for i in 1:N])
  end
  #
  # Clean up the test hdf5 file
  #
  rm("data.h5"; force=true)
end

nothing


text:
################################################################################
#   Functionality: Intakes 2 values, one from each port and outputs the sum.
#   Author: Aneta Pawelec
#   Date: 12.11.2020

function Addition_main(inPort1, inPort2, outPort1, variables)
    number1 = take!(inPort1)
    number2 = take!(inPort2)

    sum = parse(Int64 ,number1) + parse(Int64 ,number2)

    put!(outPort1, sum)
end


text:
module SimpleODEMod

using PDESolver
using SolverCommon
using ArrayViews
import ArrayViews.view
using ODLCommonTools
using SummationByParts
using PdePumiInterface
using NonlinearSolvers  # need this because not using call_nlsolver
using Jacobian
using MPI
using Utils
using Input
import ODLCommonTools: sview
export SimpleODEData, SimpleODEData_ #getMass, array3DTo1D, array1DTo3D
export evalResidual, init, run_simpleode # exported from simpleODE.jl
export ICDict              # exported from ic.jl
export ode_pre_func, ode_post_func    # exported from simpleODE_func.jl

abstract type AbstractSimpleODEData{Tsol, Tres} <: AbstractSolutionData{Tsol, Tres} end
abstract type SimpleODEData{Tsol, Tres, Tdim} <: AbstractSimpleODEData{Tsol, Tres} end

# TODO: which of these
include("types.jl")
include("check_options.jl")
include(joinpath(Pkg.dir("PDESolver"), "src/solver/debug.jl"))  # debug macro
include("simpleODE.jl")
include("common_funcs.jl")
include("ic.jl")
include("eqn_deepcopy.jl")
include("startup_func.jl")

global const PhysicsName = "SimpleODE"
register_physics(PhysicsName, SimpleODEMod, createObjects, checkOptions)

@doc """
### SimpleODEMod.calcMassMatrix

  This function calculate the mass matrix and returns it.
  Beause w are using SBP operators, the mass matrix is diagonal, so it is
  stored in a vector.

  Arguments:
    mesh: AbstractMesh
    sbp: SBP operator
    eqn: an implementation of SimpleODEData. Does not have to be fully initialized.

  Outputs:
    M: vector containing mass matrix

"""->
function calcMassMatrix(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP, eqn::SimpleODEData{Tsol, Tres, Tdim}) where {Tmsh,  Tsol, Tres, Tdim}
# calculate the (diagonal) mass matrix as a vector
# return the vector M

  M = zeros(Tmsh, mesh.numDof)
  for i=1:mesh.numEl
    for j=1:sbp.numnodes
      for k=1:mesh.numDofPerNode
        dofnum_k = mesh.dofs[k,j,i]
        # multiplication is faster than division, so do the division here
        # and then multiply solution vector times M
        M[dofnum_k] += (sbp.w[j]/mesh.jac[j,i])
      end
    end
  end

  return M

end     # end of calcMassMatrix function

# functions needed to make it compatible with the NonLinearSolvers module
function matVecA0inv(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP, eqn::SimpleODEData{Tsol, Tres, Tdim},
opts, res_arr::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tdim, Tres}

  return nothing
end

function matVecA0(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP, eqn::SimpleODEData{Tsol, Tres, Tdim}, opts,
res_arr::AbstractArray{Tsol, 3}) where {Tmsh, Tsol, Tdim, Tres}

  return nothing
end

function majorIterationCallback(itr, mesh::AbstractMesh, sbp::AbstractSBP, eqn::AbstractSimpleODEData, opts)

  return nothing
end

end # end module


text:
#printing that is not needed in Atom
import Base: show
df_show_type(io::IO, ::Type{T}) where T = dfprintln(io, crayon"red", "$T:", crayon"reset")
df_show_type(io::IO, x) = df_show_type(io, typeof(x))

function show(io::IO, block::InputData)
    s = """Block name: $(block.name)
    Block option: $(block.option)
    Block data:
    """
    dfprintln(io, s)
    dfprintln(io, string(block.data) * "\n\n")
end

show(io::IO, data::Vector{InputData}) = map(x-> show(io, x), data)

function show(io::IO, band::DFBand{T}) where T <: AbstractFloat
    df_show_type(io, band)
    string = """
    k_points of length $(length(band.k_points_cart)):
    cart:    $(band.k_points_cart[1]) -> $(band.k_points_cart[end])
    """
    !isempty(band.k_points_cryst) && (string *= "cryst:   $(band.k_points_cryst[1]) -> $(band.k_points_cryst[end])\n")
    string *= """eigvals: $(band.eigvals[1]) -> $(band.eigvals[end])
    extra:   $(band.extra)
    """
    dfprintln(io, string)
end

show(io::IO, bands::Vector{<:DFBand}) = map(x->show(io,x),bands)

function show(io::IO, job::DFJob)
    reset = crayon"reset"
    fieldns  = [:name, :local_dir, :server, :server_dir]
    fs   = string.(filter(x->!isempty(x), getfield.((job,), fieldns)))
    fns = string.(fieldns)
    lfns = maximum(length.(fns)) + maximum(length.(fs)) + 4
    line = "+"
    for i=1:div(lfns,2) + 1
        line *= "-"
    end
    dfprint(io, crayon"cyan", line[1:end-2])
    dfprint(io, crayon"cyan", "DFJOB")
    for i=1:div(lfns,2)
        line *= "-"
    end
    line *= "+"
    totlen = length(line)
    dfprintln(io, crayon"cyan", line[div(lfns, 2)+6:end],  reset)
    lfns = maximum(length.(string.(fns)))
    for (fn, f) in zip(fns, fs)
        isname = fn == "name"
        l = length(fn)
        fn *= ":"
        for i=1:lfns-l
            fn *= " "
        end
        dfprint(io, crayon"cyan", "|", reset)
        isname ? dfprint(io, " $fn ", crayon"magenta", "$f", reset) : dfprint(io, " $fn $f", reset)
        for i=1:totlen - (length(fn) + 4 + length(f))
            dfprint(io, " ")
        end
        dfprintln(io, crayon"cyan", "|", reset)
    end
    is = inputs(job)
    dfprintln(io, crayon"cyan", line, reset)
    dfprintln(io, reset,"(", crayon"green", "scheduled", reset, ", ", crayon"red", "not scheduled", reset, ")")
    ln = maximum(length.(string.(name.(is))))
    for (si, i) in enumerate(is)
        n = name(i)
        l = length(n)
        for j=1:ln-l
            n *= " "
        end
        cr = i.run ? crayon"green" : crayon"red"
        dfprint(io, cr, "\t\t$n\n")
    end
    dfprint(io, reset)
end

function show(io::IO, in::DFInput)
    df_show_type(io, in)
    s = """name  = $(in.name)
    dir   = $(in.dir)
    execs = $(join([e.exec for e in in.execs],", "))
    run   = $(in.run)
    data  = $([e.name for e in data(in)])
    $(crayon"cyan")flags$(crayon"reset"):"""
    dfprintln(io, s)
    fl = maximum(length.(string.(keys(flags(in)))))
    for (f, v) in flags(in)
        fs = string(f)
        l = length(fs)
        for i=1:fl - l
            fs *= " "
        end
        dfprint(io, crayon"cyan", "\t$fs", crayon"yellow"," => ", crayon"magenta", "$v\n")
    end
    dfprint(io, crayon"reset")
end

function show(io::IO, flag_info::QEFlagInfo{T}) where T
    df_show_type(io, flag_info)
    dfprintln(io, crayon"yellow", "name : $(flag_info.name)")
    dfprintln(io, crayon"cyan", "description:", crayon"reset")
    dfprint(io, "\t"*replace(flag_info.description, "\n" => "\n\t"))
end

function show(io::IO, flag_info::ElkFlagInfo{T}) where T
    df_show_type(io, flag_info)
    dfprintln(io, crayon"yellow", "name : $(flag_info.name)")
    dfprintln(io, crayon"green", "default:", crayon"reset")
    if flag_info.default != nothing
	    dfprintln(io, "\t" * "$(flag_info.default)")
    else
	    dfprintln(io, "\t" * "nothing")
    end
    dfprintln(io, crayon"cyan", "description:", crayon"reset")
    dfprintln(io, "\t"*replace(flag_info.description, "\n" => "\n\t"))
end

function show(io::IO, info::ElkControlBlockInfo)
    dfprintln(io, crayon"yellow", "name = $(info.name)")
    dfprintln(io, crayon"cyan", "flags:", crayon"reset")
    for flag in info.flags
	    dfprintln(io, "\t$(flag.name)")
    end
    dfprintln(io, crayon"cyan", "description:", crayon"reset")
    dfprintln(io, "\t"*replace(info.description, "\n" => "\n\t"))
end

function show(io::IO, el::Element)
	for f in fieldnames(typeof(el))
		dfprintln(io, crayon"red", "$f: ", crayon"reset", "$(getfield(el, f))")
	end
end

function show(io::IO, at::AbstractAtom{T, LT}) where {T,LT<:Length{T}}
	dfprintln(io, crayon"cyan", "$(typeof(at))")
	dfprintln(io, crayon"red", "\tname: ", crayon"reset", "$(name(at))")
	for f in fieldnames(typeof(at))[3:end-1]
		dfprintln(io, crayon"red", "\t$f: ", crayon"reset", "$(getfield(at, f))")
	end
	dfprint(io, crayon"red", "\tdftu: ", crayon"reset")
	for f in fieldnames(DFTU)
		val = getfield(dftu(at), f)
		if !isdefault(val)
			dfprint(io, "$f: $val, ")
		end
	end
	dfprintln(io)
end


text:
# Autogenerated wrapper script for LevelDB_jll for aarch64-linux-gnu-cxx11
export leveldb

JLLWrappers.@generate_wrapper_header("LevelDB")
JLLWrappers.@declare_library_product(leveldb, "libleveldb.so.1")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        leveldb,
        "lib/libleveldb.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type Sketch{T} end
abstract type FactoredSketch{T} <: Sketch{T} end

# ------------------------------------------------------------------------------
# |                               Nystrom Sketch                               |
# ------------------------------------------------------------------------------
struct NystromSketch{T} <: FactoredSketch{T}
    U::Matrix{T}
    Λ::Diagonal{T, Vector{T}}
end

#TODO: should this just wrap / generalize a NystromSketch (inc. flag for PSD)
# Specifically NystromSketch being a subtype of eigensketch
struct EigenSketch{T} <: FactoredSketch{T}
    U::Matrix{T}
    Λ::Diagonal{T, Vector{T}}
end

# Constructs Â_nys in factored form
# Â_nys = (AΩ)(ΩᵀAΩ)^†(AΩ)^ᵀ = UΛUᵀ
# [Martinsson & Tropp, Algorithm 16]
function NystromSketch(A::AbstractMatrix{T}, k::Int, r::Int; check=false) where {T <: Real}
    check && check_psd(A)
    n = size(A, 1)

    Ω = randn(n, r)
    Ω .= Array(qr(Ω).Q)
    Y = A * Ω
    ν = sqrt(n)*eps(norm(Y))                    #TODO: revisit this choice
    @. Y += ν * Ω
    B = Y / cholesky(Symmetric(Ω' * Y)).U
    U, Σ, _ = svd(B)
    Λ = Diagonal(max.(0, Σ.^2 .- ν)[1:k])

    return NystromSketch(U[:, 1:k], Λ)
end

check_input(A, ::Type{NystromSketch}) = check_psd(A)
Sketch(A, k, r, ::Type{NystromSketch}; check=false, q=0) = NystromSketch(A, k, r; check=check)

# Define basic properties
Base.eltype(::Sketch{T}) where {T} = T
Base.size(Ahat::Union{NystromSketch, EigenSketch}) = (size(Ahat.U, 1), size(Ahat.U, 1))
Base.size(Ahat::FactoredSketch, d::Int) = d <= 2 ? size(Ahat)[d] : 1
LinearAlgebra.adjoint(A::Sketch) = Adjoint(A)
LinearAlgebra.rank(Ahat::FactoredSketch) = size(Ahat.U, 2)
LinearAlgebra.svdvals(Ahat::FactoredSketch) = Ahat.Λ.diag
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, Ahat::Union{NystromSketch, EigenSketch})
    Base.showarg(stdout, Ahat, true)
    println(io, "U factor:")
    show(io, mime, Ahat.U)
    println(io, "\nEigenvalues:")
    show(io, mime, Ahat.Λ.diag)
end
#TODO: this is a hack and definitely not best practice
function Base.show(io::IO, Ahat_::Adjoint{<:Any, <:Sketch})
    print(io, "Adjoint of ")
    show(io, parent(Ahat_))
end

# Define operations for Nystrom Sketch
Matrix(Anys::Union{NystromSketch, EigenSketch}) = Anys.U*Anys.Λ*Anys.U'
LinearAlgebra.eigvals(Anys::Union{NystromSketch, EigenSketch}) = Anys.Λ.diag    # decreasing order

function LinearAlgebra.mul!(y, Anys::Union{NystromSketch, EigenSketch}, x; cache=zeros(rank(Anys)))
    length(y) != length(x) || length(y) != size(Anys, 1) && error(DimensionMismatch())
    r = rank(Anys)
    @views mul!(cache[1:r], Anys.U', x)
    @views cache[1:r] .*= Anys.Λ.diag
    @views mul!(y, Anys.U, cache[1:r])
    return nothing
end

function LinearAlgebra.:*(Anys::Union{NystromSketch, EigenSketch}, x::AbstractVector)
    n = size(Anys, 1)
    y = zeros(n)
    mul!(y, Anys, x)
    return y
end

# ------------------------------------------------------------------------------
# |                               Randomized Eigendecomposition                               |
# ------------------------------------------------------------------------------
function EigenSketch(A::AbstractMatrix{T}, k::Int, r::Int; check=false, q::Int=0) where {T <: Real}
    check && !issymmetric(A) && throw(ArgumentError("A must be symmetric"))
    k > r && throw(ArgumentError("k must be less than r"))
    Q = rangefinder(A, r; q=q)
    C = Q' * A * Q
    @. C = 0.5 * (C + C')                       #TODO: is this needed for numerics?
    Λ, V = eigen(C; sortby=x->-real(x))
    V .= real.(V)
    Λ .= real(Λ)

    # Find largest magnitude eigenspace
    p = sortperm(Λ; by=x->-abs(x))
    Λ = Λ[p][1:k]
    V = V[:, p][:,1:k]

    # Return st eigvals sorted biggest to smallest
    pp = sortperm(Λ; by=x->-x)
    Λ .= Λ[pp]
    V .= V[:, pp]

    U = Q*V
    return EigenSketch(U[:, 1:k], Diagonal(Λ))
end

check_input(A, ::Type{EigenSketch}) = issymmetric(A)
Sketch(A, k, r, ::Type{EigenSketch}; check=false, q=5) = EigenSketch(A, k, r; check=check, q=q)

# ------------------------------------------------------------------------------
# |                               Randomized SVD                               |
# ------------------------------------------------------------------------------
struct RandomizedSVD{T} <: FactoredSketch{T}
    U::Matrix{T}
    Λ::Diagonal{T, Vector{T}}
    V::Matrix{T}
end

# [Martinsson & Tropp, Algorithm 8]
# TODO: look at Tropp et al. 2019 & Halko et al. 2011 for more details on implementation
#   ‘Streaming low-rank matrix approximation with an application to scientiﬁc simulation’
#   'Finding structure with randomness: probabalistic algorithms for constructing
#    approximate matrix decompositions'
function RandomizedSVD(A::AbstractMatrix{T}, k::Int, r::Int; q::Int=0) where {T <: Real}
    k > r && throw(ArgumentError("k must be less than r"))
    Q = rangefinder(A, r; q=q)
    C = Q'*A
    Û, Σ, V = svd(C)
    U = Q*Û
    return RandomizedSVD(U[:, 1:k], Diagonal(Σ[1:k]), V[:, 1:k])
end

check_input(A, ::Type{RandomizedSVD}) = nothing
Sketch(A, k, r, ::Type{RandomizedSVD}; check=false, q=5) = RandomizedSVD(A, k, r; q=q)

# Define basic operations
Base.size(Ahat::RandomizedSVD) = (size(Ahat.U, 1), size(Ahat.V, 1))
Matrix(Ahat::RandomizedSVD) = Ahat.U*Ahat.Λ*Ahat.V'
function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, Ahat::RandomizedSVD)
    Base.showarg(stdout, Ahat, true)
    println(io, "\nU factor:")
    show(io, mime, Ahat.U)
    println(io, "\nsingular values:")
    show(io, mime, Ahat.Λ.diag)
    println(io, "\nV factor:")
    show(io, mime, Ahat.V)
end

function LinearAlgebra.mul!(y, Ahat::RandomizedSVD, x; cache=zeros(rank(Ahat)))
    length(x) != size(Ahat, 2) || length(y) != size(Ahat, 1) && error(DimensionMismatch())
    r = rank(Ahat)
    @views mul!(cache[1:r], Ahat.V', x)
    @views cache[1:r] .*= Ahat.Λ.diag
    @views mul!(y, Ahat.U, cache[1:r])
    return nothing
end

function LinearAlgebra.mul!(x, Ahat_::Adjoint{T, RandomizedSVD{T}}, y; cache=zeros(rank(Ahat))) where {T}
    Ahat = parent(Ahat_)
    length(x) != size(Ahat, 2) || length(y) != size(Ahat, 1) && error(DimensionMismatch())
    r = rank(Ahat)
    @views mul!(cache[1:r], Ahat.U', y)
    @views cache[1:r] .*= Ahat.Λ.diag
    @views mul!(x, Ahat.V, cache[1:r])
    return nothing
end

function LinearAlgebra.:*(Ahat::Union{Adjoint{T, RandomizedSVD{T}}, RandomizedSVD{T}}, x::AbstractVector) where {T}
    y = zeros(size(Ahat, 1))
    mul!(y, Ahat, x)
    return y
end

# ------------------------------------------------------------------------------
# |                             General Utilities                              |
# ------------------------------------------------------------------------------
# Power method to estimate ||A - Ahat||
function estimate_norm_E(A::AbstractMatrix{T}, Ahat::Sketch{T}; q=10, cache=nothing) where {T <: Number}
    m, n = size(A)
    if !isnothing(cache)
        u, v = cache.u, cache.v
    else
        u, v = zeros(m), zeros(n)
        cache = (Ahat_mul=zeros(min(m, n)),)
    end

    u .= randn(m)
    normalize!(u)
    v .= randn(n)
    normalize!(v)

    Ehat = Inf
    for _ in 1:q
        # u = (A - Ahat)*v
        mul!(u, Ahat, v; cache=cache.Ahat_mul)
        mul!(u, A, v, 1.0, -1.0)
        normalize!(u)

        # v = (A - Ahat)ᵀ*v
        if typeof(Ahat) <: Union{NystromSketch, EigenSketch}
            v .= u
        else
            mul!(v, Ahat', u; cache=cache.Ahat_mul)
            mul!(v, A', u, 1.0, -1.0)
            normalize!(v)
        end

        Ehat = abs(dot(u, A, v) - dot(u, Ahat, v))
    end
    return Ehat
end

# Increases rank until the approximation is sufficient
# By [Frangella et al., Prop 5.3], have that κ(P^{-1/2} * A * P^{-1/2}) ≤ (λᵣ + μ + ||E||)/μ
# TODO: Add verbose logging
# TODO: Could improve efficiency here, especially if the same sketch matrix is reused
function adaptive_sketch(
    A::AbstractMatrix{T}, r0::Int, SketchType::Type{<:Sketch};
    condition_number=false,
    ρ=1e-4,
    r_inc_factor=2.0,
    k_factor=0.9,
    tol=1e-6,
    check=false,
    q_norm=20,
    q_sketch=5,
    verbose=false
) where {T <: Real}
    check && check_input(A, SketchType)
    m, n = size(A)
    cache = (
        u=zeros(m),
        v=zeros(n),
        Ahat_mul=zeros(n)
    )
    r = r0
    Ahat = nothing
    error_metric = Inf
    while error_metric > tol && r < n
        k = round(Int, k_factor*r)
        Ahat = Sketch(A, k, r, SketchType; check=check, q=q_sketch)
        if condition_number
            error_metric = (Ahat.Λ[end] + ρ) / ρ - 1
            verbose && @info "κ = $error_metric, r = $r"
        else
            error_metric = estimate_norm_E(A, Ahat; q=q_norm, cache=cache)
            verbose && @info "||E|| = $error_metric, r = $r"
        end
        r = round(Int, r_inc_factor*r)
    end
    return Ahat
end


text:
export foundation, add_check, add_card

"""
A `foundation` consists of the four piles (one per suit) built up
starting with aces and ending with kings.
"""
struct foundation
    piles::Dict{Symbol,Vector{Card}}
    function foundation()
        p = Dict{Symbol,Vector{Card}}()
        for s in PlayingCards52.suit_list
            p[s] = Card[]
        end
        new(p)
    end
end

function length(F::foundation)::Int
    return sum(length(p[end]) for p in F.piles)
end

function string(F::foundation)::String
    result = ""
    for k = 1:4
        s = PlayingCards52.suit_list[k]
        result *= row_string(F.piles[s])
        if k < 4
            result *= "\n"
        end
    end
    return result
end

function show(io::IO, F::foundation)
    print(io, string(F))
end

"""
`add_check(F::foundation, C::Card)::Bool`
checks to see if `C` may be added to `F`.
"""
function add_check(F::foundation, PC::Card)::Bool
    s = suit(PC)
    r = rank(PC)

    pile = F.piles[s]
    if length(pile) == 0  # may only add an Ace to an empty pile
        if r == 1
            return true
        else
            return false
        end
    end

    if rank(pile[end]) + 1 == r  # may add if our rank is 1 larger than last card
        return true
    end

    return false
end

"""
`add_card(F::foundation, PC::Card)` adds the card `PC`
to `F` returning `true` if successful.
"""
function add_card(F::foundation, PC::Card)::Bool
    if !add_check(F, PC)
        return false
    end
    s = suit(PC)
    push!(F.piles[s], PC)
    return true
end

"""
`score(F::foundation)` is the sum of the squares of the lengths
of the four foundation piles.
"""
function score(F::foundation)
    sizes = length.(values(F.piles))
    return sum(sizes .^ 2)
end


text:
module RadioButtonGroups

Base.Experimental.@optlevel 1

using Redux
using CImGui
using ..RadioButtons
import ..RadioButtons: AbstractRadioButtonAction, AbstractRadioButtonState

# actions
abstract type AbstractRadioButtonGroupAction <: AbstractSyncAction end

get_label(a::AbstractRadioButtonGroupAction) = a.label

"""
    Rename(label, new_label)
Change radio button group's label to `new_label`.
Note that renaming may also change the identifier, please refer to `help?> CImGui.PushID`.
"""
struct Rename <: AbstractRadioButtonGroupAction
    label::String
    new_label::String
end

"""
    SetTriggeredTo(label, is_triggered)
Update widget's `is_triggered` state.
"""
struct SetTriggeredTo <: AbstractRadioButtonGroupAction
    label::String
    is_triggered::Bool
end

"""
    EditButtons(label::AbstractString, action)
Edit radio buttons in the radio button group.
"""
struct EditButtons{T<:AbstractRadioButtonAction} <: AbstractRadioButtonGroupAction
    label::String
    action::T
end

# state
abstract type AbstractRadioButtonGroupstate <: AbstractImmutableState end

"""
    RadioButtonGroups.State(label::AbstractString, is_active = false, is_triggered = false)
A radio button state which contains a label a.k.a the identifier, a list of radio buttons and
a flag value `is_triggered` that records the state of the latest poll events.
"""
struct State <: AbstractRadioButtonState
    label::String
    buttons::Vector{RadioButtons.State}
    is_triggered::Bool
end
State(label::AbstractString, buttons = []) = State(label, buttons, false)

get_label(s::State) = s.label
is_triggered(s::State) = s.is_triggered

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

reducer(s::State, a::Rename) = State(a.new_label, s.buttons, s.is_triggered)
reducer(s::State, a::SetTriggeredTo) = State(s.label, s.buttons, a.is_triggered)
reducer(s::State, a::EditButtons) =
    State(s.label, RadioButtons.reducer(s.buttons, a.action), s.is_triggered)

reducer(s::Dict{String,<:AbstractRadioButtonState}, a::AbstractRadioButtonGroupAction) =
    Dict(k => (get_label(v) == get_label(a) ? reducer(v, a) : v) for (k, v) in s)

# helper
"""
    RadioButtonGroup(store::AbstractStore, get_state=Redux.get_state) -> Bool
Return `true` when triggered.
`get_state` is a router function that tells how to find the target state from `store`. `chain_action` is for chaining upstream actions.
"""
function RadioButtonGroup(store::AbstractStore, get_state=Redux.get_state)
    s = get_state(store)
    is_triggered = false
    for rb in s.buttons
        flag = CImGui.RadioButton(rb.label, rb.is_active)
        dispatch!(store, EditButtons(s.label, RadioButtons.SetTriggeredTo(rb.label, flag)))
        if flag
            dispatch!(store, EditButtons(s.label, RadioButtons.OnlySetThisOneToActive(rb.label)))
            is_triggered = true
        end
        CImGui.SameLine()
    end
    CImGui.NewLine()
    is_triggered && dispatch!(store, SetTriggeredTo(s.label, is_triggered))
    return is_triggered
end

end # module