text:
struct BrowserDisplay <: Base.Multimedia.AbstractDisplay end

"""
    browser_display()
Forces JSServe.App to be displayed in a browser window that gets opened.
"""
function browser_display()
    displays = Base.Multimedia.displays
    if last(displays) isa BrowserDisplay
        @info("already in there!")
        return
    end
    # if browserdisplay is anywhere not at the last position
    # remove it!
    filter!(x-> !(x isa BrowserDisplay), displays)
    # add it to end!
    Base.pushdisplay(BrowserDisplay())
    return
end

"""
    tryrun(cmd::Cmd)

Try to run a command. Return `true` if `cmd` runs and is successful (exits with a code of `0`).
Return `false` otherwise.
"""
function tryrun(cmd::Cmd)
    try
        return success(cmd)
    catch e
        return false
    end
end

function openurl(url::String)
    if Sys.isapple()
        tryrun(`open $url`) && return
    elseif Sys.iswindows()
        tryrun(`powershell.exe start $url`) && return
    elseif Sys.isunix()
        tryrun(`xdg-open $url`) && return
        tryrun(`gnome-open $url`) && return
    end
    tryrun(`python -mwebbrowser $(url)`) && return
    # our last hope
    tryrun(`python3 -mwebbrowser $(url)`) && return
    @warn("Can't find a way to open a browser, open $(url) manually!")
end

function Base.display(::BrowserDisplay, dom::App)
    application = get_server()
    session_url = "/browser-display"
    route_was_present = route!(application, session_url) do context
        # Serve the actual content
        session = insert_session!(context.application)
        html_dom = Base.invokelatest(dom.handler, session, context.request)
        return html(page_html(session, html_dom))
    end
    # Only open url first time!
    if isempty(application.sessions)
        openurl(online_url(application, session_url))
    else
        for (id, session) in application.sessions
            evaljs(session, js"location.reload(true)")
        end
    end
    return
end


text:
##
## Join / merge
##

# Like similar, but returns a nullable array
similar_nullable{T}(dv::AbstractArray{T}, dims::@compat(Union{Int, Tuple{Vararg{Int}}})) =
    NullableArray(T, dims)

similar_nullable{T<:Nullable}(dv::AbstractArray{T}, dims::@compat(Union{Int, Tuple{Vararg{Int}}})) =
    NullableArray(eltype(T), dims)

similar_nullable{T,R}(dv::CategoricalArray{T,R}, dims::@compat(Union{Int, Tuple{Vararg{Int}}})) =
    NullableCategoricalArray(T, dims)

similar_nullable(df::AbstractDataFrame, dims::Int) =
    DataFrame(Any[similar_nullable(x, dims) for x in columns(df)], copy(index(df)))

function join_idx(left, right, max_groups)
    ## adapted from Wes McKinney's full_outer_join in pandas (file: src/join.pyx).

    # NULL group in location 0

    left_sorter, where, left_count = groupsort_indexer(left, max_groups)
    right_sorter, where, right_count = groupsort_indexer(right, max_groups)

    # First pass, determine size of result set
    tcount = 0
    rcount = 0
    lcount = 0
    for i in 1:(max_groups + 1)
        lc = left_count[i]
        rc = right_count[i]

        if rc > 0 && lc > 0
            tcount += lc * rc
        elseif rc > 0
            rcount += rc
        else
            lcount += lc
        end
    end

    # group 0 is the NULL group
    tposition = 0
    lposition = 0
    rposition = 0

    left_pos = 0
    right_pos = 0

    left_indexer = Array(Int, tcount)
    right_indexer = Array(Int, tcount)
    leftonly_indexer = Array(Int, lcount)
    rightonly_indexer = Array(Int, rcount)
    for i in 1:(max_groups + 1)
        lc = left_count[i]
        rc = right_count[i]
        if rc == 0
            for j in 1:lc
                leftonly_indexer[lposition + j] = left_pos + j
            end
            lposition += lc
        elseif lc == 0
            for j in 1:rc
                rightonly_indexer[rposition + j] = right_pos + j
            end
            rposition += rc
        else
            for j in 1:lc
                offset = tposition + (j-1) * rc
                for k in 1:rc
                    left_indexer[offset + k] = left_pos + j
                    right_indexer[offset + k] = right_pos + k
                end
            end
            tposition += lc * rc
        end
        left_pos += lc
        right_pos += rc
    end

    ## (left_sorter, left_indexer, leftonly_indexer,
    ##  right_sorter, right_indexer, rightonly_indexer)
    (left_sorter[left_indexer], left_sorter[leftonly_indexer],
     right_sorter[right_indexer], right_sorter[rightonly_indexer])
end

function sharepools{S,N}(v1::Union{CategoricalArray{S,N}, NullableCategoricalArray{S,N}},
                         v2::Union{CategoricalArray{S,N}, NullableCategoricalArray{S,N}},
                         index::Vector{S},
                         R)
    tidx1 = convert(Vector{R}, indexin(CategoricalArrays.index(v1.pool), index))
    tidx2 = convert(Vector{R}, indexin(CategoricalArrays.index(v2.pool), index))
    refs1 = zeros(R, length(v1))
    refs2 = zeros(R, length(v2))
    for i in 1:length(refs1)
        if v1.refs[i] != 0
            refs1[i] = tidx1[v1.refs[i]]
        end
    end
    for i in 1:length(refs2)
        if v2.refs[i] != 0
            refs2[i] = tidx2[v2.refs[i]]
        end
    end
    pool = CategoricalPool{S, R}(index)
    return (CategoricalArray(refs1, pool),
            CategoricalArray(refs2, pool))
end

function sharepools{S,N}(v1::Union{CategoricalArray{S,N}, NullableCategoricalArray{S,N}},
                         v2::Union{CategoricalArray{S,N}, NullableCategoricalArray{S,N}})
    index = sort(unique([levels(v1); levels(v2)]))
    sz = length(index)

    R = sz <= typemax(UInt8)  ? UInt8 :
        sz <= typemax(UInt16) ? UInt16 :
        sz <= typemax(UInt32) ? UInt32 :
                                UInt64

    # To ensure type stability during actual work
    sharepools(v1, v2, index, R)
end

sharepools{S,N}(v1::Union{CategoricalArray{S,N}, NullableCategoricalArray{S,N}},
                v2::AbstractArray{S,N}) =
    sharepools(v1, oftype(v1, v2))

sharepools{S,N}(v1::AbstractArray{S,N},
                v2::Union{CategoricalArray{S,N}, NullableCategoricalArray{S,N}}) =
    sharepools(oftype(v2, v1), v2)

# TODO: write an optimized version for (Nullable)CategoricalArray
function sharepools(v1::AbstractArray,
                    v2::AbstractArray)
    ## Return two categorical arrays that share the same pool.

    ## TODO: allow specification of R
    R = CategoricalArrays.DefaultRefType
    refs1 = Array(R, size(v1))
    refs2 = Array(R, size(v2))
    poolref = Dict{promote_type(eltype(v1), eltype(v2)), R}()
    maxref = 0

    # loop through once to fill the poolref dict
    for i = 1:length(v1)
        if !_isnull(v1[i])
            poolref[v1[i]] = 0
        end
    end
    for i = 1:length(v2)
        if !_isnull(v2[i])
            poolref[v2[i]] = 0
        end
    end

    # fill positions in poolref
    pool = sort(collect(keys(poolref)))
    i = 1
    for p in pool
        poolref[p] = i
        i += 1
    end

    # fill in newrefs
    zeroval = zero(R)
    for i = 1:length(v1)
        if _isnull(v1[i])
            refs1[i] = zeroval
        else
            refs1[i] = poolref[v1[i]]
        end
    end
    for i = 1:length(v2)
        if _isnull(v2[i])
            refs2[i] = zeroval
        else
            refs2[i] = poolref[v2[i]]
        end
    end

    pool = CategoricalPool(pool)
    return (NullableCategoricalArray(refs1, pool),
            NullableCategoricalArray(refs2, pool))
end

function sharepools(df1::AbstractDataFrame, df2::AbstractDataFrame)
    # This method exists to allow merge to work with multiple columns.
    # It takes the columns of each DataFrame and returns a categorical array
    # with a merged pool that "keys" the combination of column values.
    # The pools of the result don't really mean anything.
    dv1, dv2 = sharepools(df1[1], df2[1])
    # use UInt32 instead of the minimum integer size chosen by sharepools
    # since the number of levels can be high
    refs1 = Vector{UInt32}(dv1.refs)
    refs2 = Vector{UInt32}(dv2.refs)
    # the + 1 handles nulls
    refs1[:] += 1
    refs2[:] += 1
    ngroups = length(levels(dv1)) + 1
    for j = 2:ncol(df1)
        dv1, dv2 = sharepools(df1[j], df2[j])
        for i = 1:length(refs1)
            refs1[i] += (dv1.refs[i]) * ngroups
        end
        for i = 1:length(refs2)
            refs2[i] += (dv2.refs[i]) * ngroups
        end
        ngroups *= length(levels(dv1)) + 1
    end
    # recode refs1 and refs2 to drop the unused column combinations and
    # limit the pool size
    sharepools(refs1, refs2)
end

"""
Join two DataFrames

```julia
join(df1::AbstractDataFrame,
     df2::AbstractDataFrame;
     on::Union{Symbol, Vector{Symbol}} = Symbol[],
     kind::Symbol = :inner)
```

### Arguments

* `df1`, `df2` : the two AbstractDataFrames to be joined

### Keyword Arguments

* `on` : a Symbol or Vector{Symbol}, the column(s) used as keys when
  joining; required argument except for `kind = :cross`

* `kind` : the type of join, options include:

  - `:inner` : only include rows with keys that match in both `df1`
    and `df2`, the default
  - `:outer` : include all rows from `df1` and `df2`
  - `:left` : include all rows from `df1`
  - `:right` : include all rows from `df2`
  - `:semi` : return rows of `df1` that match with the keys in `df2`
  - `:anti` : return rows of `df1` that do not match with the keys in `df2`
  - `:cross` : a full Cartesian product of the key combinations; every
    row of `df1` is matched with every row of `df2`

Null values are filled in where needed to complete joins.

### Result

* `::DataFrame` : the joined DataFrame

### Examples

```julia
name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])

join(name, job, on = :ID)
join(name, job, on = :ID, kind = :outer)
join(name, job, on = :ID, kind = :left)
join(name, job, on = :ID, kind = :right)
join(name, job, on = :ID, kind = :semi)
join(name, job, on = :ID, kind = :anti)
join(name, job, kind = :cross)
```

"""
function Base.join(df1::AbstractDataFrame,
                   df2::AbstractDataFrame;
                   on::@compat(Union{Symbol, Vector{Symbol}}) = Symbol[],
                   kind::Symbol = :inner)
    if kind == :cross
        if on != Symbol[]
            throw(ArgumentError("Cross joins don't use argument 'on'."))
        end
        return crossjoin(df1, df2)
    elseif on == Symbol[]
        throw(ArgumentError("Missing join argument 'on'."))
    end

    dv1, dv2 = sharepools(df1[on], df2[on])

    left_idx, leftonly_idx, right_idx, rightonly_idx =
        join_idx(dv1.refs, dv2.refs, length(dv1.pool))

    if kind == :inner
        df2w = without(df2, on)

        left = df1[left_idx, :]
        right = df2w[right_idx, :]

        return hcat!(left, right)
    elseif kind == :left
        df2w = without(df2, on)

        left = df1[[left_idx; leftonly_idx], :]
        right = vcat(df2w[right_idx, :],
                     similar_nullable(df2w, length(leftonly_idx)))

        return hcat!(left, right)
    elseif kind == :right
        df1w = without(df1, on)

        left = vcat(df1w[left_idx, :],
                    similar_nullable(df1w, length(rightonly_idx)))
        right = df2[[right_idx; rightonly_idx], :]

        return hcat!(left, right)
    elseif kind == :outer
        df1w, df2w = without(df1, on),  without(df2, on)

        mixed = hcat!(df1[left_idx, :], df2w[right_idx, :])
        leftonly = hcat!(df1[leftonly_idx, :],
                         similar_nullable(df2w, length(leftonly_idx)))
        rightonly = hcat!(similar_nullable(df1w, length(rightonly_idx)),
                          df2[rightonly_idx, :])

        return vcat(mixed, leftonly, rightonly)
    elseif kind == :semi
        df1[unique(left_idx), :]
    elseif kind == :anti
        df1[leftonly_idx, :]
    else
        throw(ArgumentError("Unknown kind of join requested"))
    end
end

function crossjoin(df1::AbstractDataFrame, df2::AbstractDataFrame)
    r1, r2 = size(df1, 1), size(df2, 1)
    cols = Any[[Compat.repeat(c, inner=r2) for c in columns(df1)];
            [Compat.repeat(c, outer=r1) for c in columns(df2)]]
    colindex = merge(index(df1), index(df2))
    DataFrame(cols, colindex)
end

ext: jl


text:
# struct containing api key and urls
struct wmata_API
    api_key::String
    station_list_url::String
    station_timings_url::String
    rail_predictions_url::String
    paths_url::String
    station_to_station_url::String
    train_positions_url::String
    rail_incidents_url::String
end


text:
using Base.Test
import ArchGDAL; const AG = ArchGDAL

@testset "Testing StyleTable Methods" begin
    AG.createstylemanager() do sm
        @test AG.initialize!(sm) == true
        @test AG.npart(sm) == 0
        @test AG.initialize!(sm, "PEN(w:2px,c:#000000,id:\"mapinfo-pen-2,ogr-pen-0\")") == true
        AG.getpart(sm, 0) do st
            @test AG.getstylestring(st) == "PEN(w:2px,c:#000000,id:\"mapinfo-pen-2,ogr-pen-0\")"
            @test AG.getrgba(st, "#123456") == (18, 52, 86, 255)
        end
        @test AG.npart(sm) == 1
        @test AG.addstyle!(sm, "name1", "style1") == false
        @test AG.npart(sm) == 1
        @test AG.addstyle!(sm, "name2") == false
        @test AG.npart(sm) == 1

        AG.createstyletool(GDAL.OGRSTCBrush) do st
            @test AG.gettype(st) == GDAL.OGRSTCBrush
            @test AG.getunit(st) == GDAL.OGRSTUMM
            AG.setunit!(st, GDAL.OGRSTUPixel, 2.0)
            @test AG.getunit(st) == GDAL.OGRSTUPixel

            AG.setparam!(st, 0, 0)
                @test AG.asint(st, 0) == 0
                @test AG.asstring(st, 0) == "0"
                @test AG.asdouble(st, 0) == 0
            AG.setparam!(st, 1, 12)
                @test AG.asint(st, 1) == 12
                @test AG.asstring(st, 1) == "12"
                @test AG.asdouble(st, 1) == 12
            AG.setparam!(st, 2, "foo")
                @test AG.asstring(st, 2) == "foo"
            AG.setparam!(st, 3, 0.5)
                @test AG.asdouble(st, 3) ≈ 0.5

            @test AG.npart(sm) == 1
            AG.addpart!(sm, st)
            @test AG.npart(sm) == 2
            @test AG.npart(sm, "some stylestring") == 1
        end

        AG.createstyletable() do stbl
            AG.addstyle!(stbl, "name1", "style1")
            AG.addstyle!(stbl, "name2", "style2")
            AG.addstyle!(stbl, "name3", "style3")
            AG.addstyle!(stbl, "name4", "style4")
            @test AG.find(stbl, "name3") == "style3"
            @test AG.laststyle(stbl) == ""
            @test AG.nextstyle(stbl) == "style1"
            @test AG.nextstyle(stbl) == "style2"
            AG.resetreading!(stbl)
            @test AG.nextstyle(stbl) == "style1"
            AG.savestyletable(stbl, "tmp/styletable.txt")
        end
        AG.createstyletable() do stbl
            AG.loadstyletable!(stbl, "tmp/styletable.txt")
            @test AG.find(stbl, "name3") == "style3"
            @test AG.laststyle(stbl) == ""
            @test AG.nextstyle(stbl) == "style1"
            @test AG.nextstyle(stbl) == "style2"
            AG.resetreading!(stbl)
            @test AG.nextstyle(stbl) == "style1"
        end
        rm("tmp/styletable.txt")
    end
end

# Untested: initialize!(stylemanager::StyleManager, feature::Feature)


text:
using DustExtinction
import DustExtinction: redden, deredden

export redden,
       redden!,
       deredden,
       deredden!

"""
    redden!(::AbstractSpectrum, Av; Rv = 3.1, law = ccm89)

In-place version of [`redden`](@ref)
"""
function redden!(spec::T, Av; Rv = 3.1, law = ccm89) where {T <: AbstractSpectrum}
    @. spec.flux = redden(spec.flux, spec.wave, Av, Rv = Rv, law = law)
    return spec
end

"""
    redden(::AbstractSpectrum, Av; Rv = 3.1, law = ccm89)

Redden a spectrum using common color laws provided by [DustExtinction.jl](https://github.com/juliaastro/dustextinction.jl). `Av` is the total extinction, `Rv` is the selective extinction (3.1 is a common value for the Milky Way) and `law` is the color law to use for determining the extinction.
"""
function redden(spec::AbstractSpectrum, Av; Rv = 3.1, law = ccm89)
    tmp_spec = deepcopy(spec)
    redden!(tmp_spec, Av, Rv = Rv, law = law)
    return tmp_spec
end

"""
    deredden!(::AbstractSpectrum, Av; Rv = 3.1, law = ccm89)

In-place version of [`deredden`](@ref)
"""
function deredden!(spec::AbstractSpectrum, Av; Rv = 3.1, law = ccm89)
    @. spec.flux = deredden(spec.flux, spec.wave, Av, Rv = Rv, law = law)
    return spec
end

"""
    deredden(::AbstractSpectrum, Av; Rv = 3.1, law = ccm89)

Deredden a spectrum using common color laws provided by [DustExtinction.jl](https://github.com/juliaastro/dustextinction.jl). `Av` is the total extinction, `Rv` is the selective extinction (3.1 is a common value for the Milky Way) and `law` is the color law to use for determining the extinction.
"""
function deredden(spec::AbstractSpectrum, Av; Rv = 3.1, law = ccm89)
    tmp_spec = deepcopy(spec)
    deredden!(tmp_spec, Av, Rv = Rv, law = law)
    return tmp_spec
end


text:
using PackageCompiler

create_sysimage([
		:ArgParse,
		:BSON,
		:CSV,
		:DSP,
		:DataFrames,
		:FFTW,
		:Flux,
		:MP3,
		:ProgressMeter,
		:UnicodePlots,
	];
	project=joinpath(@__DIR__, ".."),
	sysimage_path=joinpath(@__DIR__, "../sysimage.so"))

max_stars_repo_name: gabriel-ss/music-autotagger

max_issues_repo_name: gabriel-ss/music-autotagger

max_forks_repo_name: gabriel-ss/music-autotagger


text:
export Ac_mul_B!

function Ac_mul_B!( alpha::Float64,
                    A::SparseMatrixCSC{Float64,Int},
                    x::Array{Float64},
                    beta::Float64,
                    y::Array{Float64},
                    nthreads::Int64=0 )
# Real:  y = beta*y  +  alpha * A'*x

   if nthreads == 0
      #Ac_mul_B!( alpha, A, x, beta, y )
	  mul!(y,adjoint(A),x, alpha, beta)
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

	n,m  = size(A)
   nvec = size(x,2)

   if size(x,1) != n || size(y,1) != m
      throw(DimensionMismatch("size(x) != n || size(y) != m"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:ac_mul_b_rr_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Float64}),
                Ref(nthreads), Ref(nvec), Ref(m), Ref(n),    Ref(alpha),   Ref(beta),              A.nzval,      A.rowval,   A.colptr,   x,   y);

end  # function Ac_mul_B!

#------------------------------------------------------------------------------

function Ac_mul_B!( alpha::ComplexF64,
                    A::SparseMatrixCSC{Float64,Int},
                    x::Array{ComplexF64},
                    beta::ComplexF64,
                    y::Array{ComplexF64},
                    nthreads::Int64=0 )
# Real, Complex A:  y = beta*y  +  alpha * A'*x

   if nthreads == 0
      mul!(y,adjoint(A),x, alpha, beta)
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

	n,m  = size(A)
   nvec = size(x,2)

   if size(x,1) != n || size(y,1) != m
      throw(DimensionMismatch("length(x) != n || length(y) != m"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:ac_mul_b_rc_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{ComplexF64}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{ComplexF64}),
                   Ref(nthreads), Ref(nvec), Ref(m), Ref(n),     Ref(alpha),   Ref(beta),              A.nzval,      A.rowval,   A.colptr,   convert(Ptr{ComplexF64}, pointer(x)),  convert(Ptr{ComplexF64}, pointer(y)));

end  # function Ac_mul_B!

#------------------------------------------------------------------------------

function Ac_mul_B!( alpha::ComplexF64,
                    A::SparseMatrixCSC{ComplexF64,Int},
                    x::Array{ComplexF64},
                    beta::ComplexF64,
                    y::Array{ComplexF64},
                    nthreads::Int64=0 )
# Complex:  y = beta*y  +  alpha * A'*x

   if nthreads == 0
      mul!(y,adjoint(A),x, alpha, beta) #Base.Ac_mul_B!( alpha, A, x, beta, y )
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

   n,m  = size(A)
   nvec = size(x,2)

   if size(x,1) != n || size(y,1) != m
      throw(DimensionMismatch("length(x) != n || length(y) != m"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:ac_mul_b_cc_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{ComplexF64}, Ptr{ComplexF64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{ComplexF64}),
                   Ref(nthreads), Ref(nvec), Ref(m), Ref(n),     Ref(alpha),   Ref(beta),              convert(Ptr{ComplexF64}, pointer(A.nzval)),      A.rowval,   A.colptr,   convert(Ptr{ComplexF64}, pointer(x)),  convert(Ptr{ComplexF64}, pointer(y)));

end  # function Ac_mul_B!

function Ac_mul_B!( alpha::ComplexF32,
                    A::SparseMatrixCSC{ComplexF32,Int64},
                    x::Array{ComplexF32},
                    beta::ComplexF32,
                    y::Array{ComplexF32},
                    nthreads::Int64=0 )
# Complex:  y = beta*y  +  alpha * A'*x

   if nthreads == 0
      mul!(y,adjoint(A),x, alpha, beta) #Base.Ac_mul_B!( alpha, A, x, beta, y )
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

   n,m  = size(A)
   nvec = size(x,2)

   if size(x,1) != n || size(y,1) != m
      throw(DimensionMismatch("length(x) != n || length(y) != m"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:ac_mul_b_cc_short_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF32}, Ptr{ComplexF32}, Ptr{ComplexF32}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF32}, Ptr{ComplexF32}),
                   Ref(nthreads), Ref(nvec), Ref(m), Ref(n),     Ref(alpha),   Ref(beta),              convert(Ptr{ComplexF32}, pointer(A.nzval)),      A.rowval,   A.colptr,   convert(Ptr{ComplexF32}, pointer(x)),  convert(Ptr{ComplexF32}, pointer(y)));

end  # function Ac_mul_B!

function Ac_mul_B!( alpha::ComplexF32,
                    A::SparseMatrixCSC{Float32,Int64},
                    x::Array{ComplexF32},
                    beta::ComplexF32,
                    y::Array{ComplexF32},
                    nthreads::Int64=0 )
# Complex:  y = beta*y  +  alpha * A'*x

   if nthreads == 0
      mul!(y,adjoint(A),x, alpha, beta) #Base.Ac_mul_B!( alpha, A, x, beta, y )
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

   n,m  = size(A)
   nvec = size(x,2)

   if size(x,1) != n || size(y,1) != m
      throw(DimensionMismatch("length(x) != n || length(y) != m"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:ac_mul_b_rc_short_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF32}, Ptr{ComplexF32}, Ptr{Float32}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF32}, Ptr{ComplexF32}),
                   Ref(nthreads), Ref(nvec), Ref(m), Ref(n),     Ref(alpha),   Ref(beta),  A.nzval,      A.rowval,   A.colptr,   convert(Ptr{ComplexF32}, pointer(x)),  convert(Ptr{ComplexF32}, pointer(y)));

end  # function Ac_mul_B!

function Ac_mul_B!( alpha::ComplexF64,
                    A::SparseMatrixCSC{ComplexF32,Int64},
                    x::Array{ComplexF64},
                    beta::ComplexF64,
                    y::Array{ComplexF64},
                    nthreads::Int64=0 )
# Complex:  y = beta*y  +  alpha * A'*x

   if nthreads == 0
      mul!(y,adjoint(A),x, alpha, beta) #Base.Ac_mul_B!( alpha, A, x, beta, y )
      return
   elseif nthreads < 1
      throw(ArgumentError("nthreads < 1"))
   end

   n,m  = size(A)
   nvec = size(x,2)

   if size(x,1) != n || size(y,1) != m
      throw(DimensionMismatch("length(x) != n || length(y) != m"))
   elseif size(y,2) != nvec
      throw(DimensionMismatch("length(y,2) != nvec"))
   end

	p  = ccall( (:ac_mul_b_cc_mixed_, spmatveclib),
		 Int64, ( Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{ComplexF64}, Ptr{ComplexF32}, Ptr{Int64}, Ptr{Int64}, Ptr{ComplexF64}, Ptr{ComplexF64}),
                   Ref(nthreads), Ref(nvec), Ref(m), Ref(n),     Ref(alpha),   Ref(beta),              convert(Ptr{ComplexF32}, pointer(A.nzval)),      A.rowval,   A.colptr,   convert(Ptr{ComplexF64}, pointer(x)),  convert(Ptr{ComplexF64}, pointer(y)));

end  # function Ac_mul_B!


text:
using ParallelAnalysis
using Test, Random

@testset "ParallelAnalysis.jl" begin
    include("Test1_polychoric.jl")
    include("Test2_parallel.jl")
end


text:
using cfgrib
using DataStructures

#  Dummy offsets for testing
function dummy_offsets()
    message1 = (
        int    = 10,
        array  = [1,2,3],
        string = "potato",
        ignored = "☹",
        absent = missing,
        one_here = 20
    )

    message2 = (
        int    = 10,
        array  = [1,2,3,5],
        string = "cabbage",
        absent = missing,
        one_here = missing
    )

    dummy_offsets = OrderedDict(
        message1 => 0,
        message2 => 10_000
    )

    return dummy_offsets
end

#  Dummy FileIndex for testing
function dummy_file_index()
    fileindex = cfgrib.FileIndex()

    fileindex.grib_path = "./dummy-path.grib"

    fileindex.index_keys = ["int", "array", "string", "absent", "one_here"]
    fileindex.offsets = collect(pairs(dummy_offsets()))

    return fileindex
end

dfi = dummy_file_index()

#  When new fields are added, add any new required tests and then
#  add the field to this list
@test fieldnames(cfgrib.FileIndex) == (
    :allowed_protocol_version, :grib_path, :index_path, :index_keys, :offsets,
    :message_lengths, :header_values, :filter_by_keys
)

#  If this fails, check that any changes to the code are
#  reflected in the tests
@test dfi.allowed_protocol_version == v"0.0.0"

#  These should not be initialised yet
@test ! isdefined(dfi, :index_path)
@test ! isdefined(dfi, :header_values)
@test ! isdefined(dfi, :filter_by_keys)

#  Will fail on version change
cfgrib.index_path!(dfi)
@test dfi.index_path == "../dummy-path.grib.60099cfb35e25e30.idx"

cfgrib.get_header_values!(dfi)
expected_header_values = OrderedDict(
    "int"      => [10],
    "array"    => Array{Int64,1}[[1, 2, 3], [1, 2, 3, 5]],
    "string"   => ["potato", "cabbage"],
    "absent"   => Missing[missing],
    "one_here" => Union{Missing, Int64}[20, missing]
)
@test isequal(dfi.header_values, expected_header_values)

@testset begin
    test_file = joinpath(dir_testfiles, "era5-levels-members.grib")

    index = cfgrib.FileIndex(
        test_file,
        cfgrib.ALL_KEYS
    )

    message = cfgrib.first(index)

    cfgrib.filter!(index, paramId=130)
end

#  Not implemented yet
@test_skip cfgrib.save_indexfile!(dfi)
@test_skip cfgrib.from_indexfile!(dfi)


text:
function add_additional_state_provider(obj::FieldAbstractPropagator, arg0::FieldAdditionalStateProvider)
    return jcall(obj, "addAdditionalStateProvider", void, (FieldAdditionalStateProvider,), arg0)
end

function add_event_detector(obj::FieldAbstractAnalyticalPropagator, arg0::FieldEventDetector)
    return jcall(obj, "addEventDetector", void, (FieldEventDetector,), arg0)
end

function clear_events_detectors(obj::FieldAbstractAnalyticalPropagator)
    return jcall(obj, "clearEventsDetectors", void, ())
end

function get_additional_state_providers(obj::FieldAbstractPropagator)
    return jcall(obj, "getAdditionalStateProviders", List, ())
end

function get_attitude_provider(obj::FieldAbstractPropagator)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_events_detectors(obj::FieldAbstractAnalyticalPropagator)
    return jcall(obj, "getEventsDetectors", Collection, ())
end

function get_field(obj::FieldAbstractPropagator)
    return jcall(obj, "getField", Field, ())
end

function get_frame(obj::FieldAbstractPropagator)
    return jcall(obj, "getFrame", Frame, ())
end

function get_generated_ephemeris(obj::FieldAbstractAnalyticalPropagator)
    return jcall(obj, "getGeneratedEphemeris", FieldBoundedPropagator, ())
end

function get_initial_state(obj::FieldAbstractPropagator)
    return jcall(obj, "getInitialState", FieldSpacecraftState, ())
end

function get_managed_additional_states(obj::FieldAbstractPropagator)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_mode(obj::FieldAbstractPropagator)
    return jcall(obj, "getMode", jint, ())
end

function get_pv_coordinates(obj::FieldAbstractPropagator, arg0::FieldAbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (FieldAbsoluteDate, Frame), arg0, arg1)
end

function get_pv_provider(obj::FieldAbstractAnalyticalPropagator)
    return jcall(obj, "getPvProvider", FieldPVCoordinatesProvider, ())
end

function is_additional_state_managed(obj::FieldAbstractPropagator, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function propagate(obj::FieldAbstractAnalyticalPropagator, arg0::FieldAbsoluteDate, arg1::FieldAbsoluteDate)
    return jcall(obj, "propagate", FieldSpacecraftState, (FieldAbsoluteDate, FieldAbsoluteDate), arg0, arg1)
end

function propagate(obj::FieldAbstractPropagator, arg0::FieldAbsoluteDate)
    return jcall(obj, "propagate", FieldSpacecraftState, (FieldAbsoluteDate,), arg0)
end

function reset_initial_state(obj::FieldAbstractPropagator, arg0::FieldSpacecraftState)
    return jcall(obj, "resetInitialState", void, (FieldSpacecraftState,), arg0)
end

function set_attitude_provider(obj::FieldAbstractPropagator, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function set_ephemeris_mode(obj::FieldAbstractPropagator)
    return jcall(obj, "setEphemerisMode", void, ())
end

function set_master_mode(obj::FieldAbstractPropagator, arg0::FieldOrekitStepHandler)
    return jcall(obj, "setMasterMode", void, (FieldOrekitStepHandler,), arg0)
end

function set_master_mode(obj::FieldAbstractPropagator, arg0::RealFieldElement, arg1::FieldOrekitFixedStepHandler)
    return jcall(obj, "setMasterMode", void, (RealFieldElement, FieldOrekitFixedStepHandler), arg0, arg1)
end

function set_slave_mode(obj::FieldAbstractPropagator)
    return jcall(obj, "setSlaveMode", void, ())
end


text:
using Unitful, ImagineInterface, Imagine, Statistics
using Test

srate = 5000*inv(Unitful.s) #this is the max AO rate for the device
dummy = rigtemplate("dummy-6002"; sample_rate = srate)
pos = getpositioners(dummy)[1]
ao2 = getoutputs(getanalog(dummy))[2]
#pos_mon = ImagineInterface.getpositionermonitors(dummy)[1]
ai = getinputs(getanalog(dummy))[3] #AI2. note pos_mon currently has a problem with negative voltages
ai2 = getinputs(getanalog(dummy))[4] #AI3

#error()

nsamps_on = 10001
nsamps_off = 10001
nsamps_back_on = 10001
nsamps = nsamps_on+nsamps_off+nsamps_back_on
append!(pos, "on", fill(3.0*Unitful.V, nsamps_on))
append!(pos, "off", fill(0.0*Unitful.V, nsamps_off))
append!(pos, "back_on", fill(3.0*Unitful.V, nsamps_back_on))
#append!(ao2, "on2", fill(3.0*Unitful.V, nsamps_on))
append!(ao2, "off2", fill(0.0*Unitful.V, nsamps_off))
append!(ao2, "off2")
append!(ao2, "off2")
#append!(ao2, "back_on2", fill(3.0*Unitful.V, nsamps_back_on))
base_name = "" #"test_imaginectrl5"
aos = [pos;ao2]
ais = [ai;ai2]
#aos = [pos]
#ais = [ai]

#hardware setup:
# PFI1 <-> P0.0
# PFI2 <-> P0.0
# AO0 <-> AI2

sigs = run_imagine("",
		   vcat(aos,ais);
		   ai_trig_dest = "PFI1",
		   ao_trig_dest = "PFI0",
		   trigger_source = "Port0/Line0",
		   sync_clocks = false,
		   run_locally = false,
		   skip_validation = true)

@test Statistics.cor(ustrip.(get_samples(sigs[1]; sampmap=:volts)), ustrip.(get_samples(pos, sampmap=:volts))) >= 0.99

#using Plots
#plot(ustrip(get_samples(sigs[1])))

##for testing
#trigger_term_ao = "PFI0"
#trigger_term_ai = "PFI1"
#Imagine.ImagineWorker._set_device("Dev1/")
#tsk = Imagine.ImagineWorker.prepare_ao(aos, 1000, trigger_term_ao, "")
#tsk2 = Imagine.ImagineWorker.prepare_ai(ais, 1000, 3003, trigger_term_ai; clock_source="")
##After much pain, realized that can only synchronize start of AO and AI acq on usb 6xxx devices by using two different pfi channels
##NIDAQ.catch_error(NIDAQ.ExportSignal(tsk.th, NIDAQ.Val_StartTrigger, b"/Dev1/ai/StartTrigger"))
#rr2 = record_signals(base_name, ais, nsamps; trigger_terminal = trigger_term_ai, run_locally = false)
#rr = output_signals(aos; trigger_terminal = trigger_term_ao, run_locally = false)
#sleep(10.0)
#print("delivering pulses\n")
#Imagine.ttl_pulse(;line_name="Port0/line0") #P0.0 is wired to PFI0
#(proc_id2, ai_sigs) = fetch(rr2)
#ai_recs = fetch(ai_sigs[1])
#proc_id1 = fetch(rr[2][1])
#
#print("Freeing workers...\n")
#free_workers(proc_id1)
#free_workers(proc_id2)
##print("Loading ai file...\n")
##aidata = ImagineInterface.parse_ai(base_name*".ai", map(daq_channel, ais), "dummy-6002", srate)
#using Plots
#plot(ustrip(get_samples(sigs[1])))


text:
# AUTO GENERATED FILE - DO NOT EDIT

export vtk_imagedata

"""
    vtk_imagedata(;kwargs...)
    vtk_imagedata(children::Any;kwargs...)
    vtk_imagedata(children_maker::Function;kwargs...)

An ImageData component.
ImageData is exposing a vtkImageData to a downstream filter
It takes the following set of properties:
  - dimensions: [nx, ny, nz],
  - origin: [0, 0, 0]
  - spacing: [1, 1, 1]
  - direction: [
      1, 0, 0,
      0, 1, 0,
      0, 0, 1
    ]
Keyword arguments:
- `children` (Array of a list of or a singular dash component, string or numbers | a list of or a singular dash component, string or number; optional)
- `id` (String; optional): The ID used to identify this component.
- `port` (Real; optional): downstream connection port
- `dimensions` (Array of Reals; optional): Number of points along x, y, z
- `spacing` (Array of Reals; optional): Spacing along x, y, z between points in world coordinates
- `origin` (Array of Reals; optional): World coordinate of the lower left corner of your vtkImageData (i=0, j=0, k=0).
- `direction` (Array of Reals; optional): 3x3 matrix use to orient the image data
"""
function vtk_imagedata(; kwargs...)
        available_props = Symbol[:children, :id, :port, :dimensions, :spacing, :origin, :direction]
        wild_props = Symbol[]
        return Component("vtk_imagedata", "ImageData", "dash_vtk", available_props, wild_props; kwargs...)
end

vtk_imagedata(children::Any; kwargs...) = vtk_imagedata(;kwargs..., children = children)
vtk_imagedata(children_maker::Function; kwargs...) = vtk_imagedata(children_maker(); kwargs...)

max_stars_repo_name: ericbdaniels/dash-vtk

max_issues_repo_name: ericbdaniels/dash-vtk

max_forks_repo_name: ericbdaniels/dash-vtk


text:
#tracywidom.m
#Algorithm 1.3 of Random Eigenvalues by Alan Edelman

#Theory:      Compute and plot the Tracy-Widom distribution
#See also:    RandomMatrices.TracyWidom(t) for tn=-8

##Parameters
t0=5.         # right endpoint
tn=-8.        # left  endpoint

##Theory: The differential equation solver
#Compute Tracy-Widom distribution
using ODE
include("gradient.jl")
function tracywidom(t0::Float64, tn::Float64)
    function deq(t::Float64, y::Vector{Float64})
        yout = [y[2]; t*y[1]+2y[1]^3; y[4]; y[1]^2]
    end

    y0=[airy(t0); airy(1, t0); 0; airy(t0)^2] # Initial conditions
    (t, y)=ode23(deq, [t0, tn], y0)           # Solve the ODE
    F2=exp(-y[:,3])                           # the distribution
    f2=gradient(F2, t)                        # the density
    return (t, f2)
end
(t, f2) = tracywidom(t0, tn)

## Plot
using Winston
p = FramedPlot()
add(p, Curve(t, f2, "linewidth", 2))
if isinteractive()
    Winston.display(p)
else
    file(p, "tracywidom.png")
end


text:
#create all required objects
smhpolicy1 = policy(:mh,1)
smhpolicy2 = policy(:mh,1,store=:all)
smhrunner1 = runner(smhpolicy1,rniter1;numburnin = rnburnin1)
smhrunner2 = runner(smhpolicy2,rniter1,1;numburnin = rnburnin1)
smhindicator1,smhtunerstate1,smhchain1 = GeneralizedMetropolisHastings.createcommon(smhrunner1,rtuner1,rnparas1,smhnprops1,smhnprops1)
smhindicator2,smhtunerstate2,smhchain2 = GeneralizedMetropolisHastings.createcommon(smhrunner1,rtuner2,rnparas1,smhnprops1,smhnprops1)
smhsamplerstate1 = samplerstate(rsampler1,smhnprops1,Float64,Float64)
smhsamplerstate2 = samplerstate(rsampler1,smhnprops1,Float64,Float64)
@test numproposals(smhindicator1) == 1 && numsamples(smhindicator1) == 1
@test indicatorsamples(smhindicator1) == [0,0]
@test indicatorsamples(smhindicator2) == [0,0]
@test numtunesteps(smhtunerstate1) == 2
@test numparas(smhchain1) == rnparas1 && numsamples(smhchain1) == 4

#test the preparenext! function
sv1 = [1.0]
sll1 = [-100.0]
sv2 = [3.0]
sll2 = [-50.0]
GMHRunnersTest.setfrom(smhsamplerstate1,sv1,sll1) ; GMHRunnersTest.setproposals(smhsamplerstate1,sv2,sll2) ; GMHRunnersTest.setindicator(smhindicator1,[1.0,0.0],[2,1])
preparenext!(smhrunner1,rmodel1,smhsamplerstate1,smhindicator1)
GMHRunnersTest.testfrom(smhsamplerstate1,sv2,sll2) #in this condition, value should have been copied
@test_approx_eq mean(smhsamplerstate1.density.distribution) sv2 #and the normal distribution of the mh sampler should also have been conditioned
GMHRunnersTest.setfrom(smhsamplerstate1,sv1,sll1) ; GMHRunnersTest.setproposals(smhsamplerstate1,sv2,sll2) ; GMHRunnersTest.setindicator(smhindicator1,[0.0,1.0],[2,2])
preparenext!(smhrunner1,rmodel1,smhsamplerstate1,smhindicator1)
GMHRunnersTest.testfrom(smhsamplerstate1,sv1,sll1) #in this condition, value should have stayed the same
@test_approx_eq mean(smhsamplerstate1.density.distribution) sv2 #and the normal distribution should have been left unchanged from previous test

#test storing a value into the chain
GMHRunnersTest.setfrom(smhsamplerstate1,sv1,sll1) ; GMHRunnersTest.setproposals(smhsamplerstate1,sv2,sll2) ; GMHRunnersTest.setindicator(smhindicator1,[1.0,0.0],[2,1])
store!(smhrunner1,smhsamplerstate1,smhindicator1,smhchain1)
GMHRunnersTest.teststore(smhchain1,sv2,sll2[1],1,1,1)
GMHRunnersTest.setfrom(smhsamplerstate1,sv1,sll1) ; GMHRunnersTest.setproposals(smhsamplerstate1,sv2,sll2) ; GMHRunnersTest.setindicator(smhindicator1,[0.0,1.0],[2,2])
store!(smhrunner1,smhsamplerstate1,smhindicator1,smhchain1)
GMHRunnersTest.teststore(smhchain1,sv1,sll1[1],2,1,2)
@test smhchain1.values == [sv2[1] sv1[1] 0.0 0.0] #test the exact values

#iterate once
println("============================")
println("Testing of iterate! function")
println("============================")
#reset the samplerstates used below
smhsamplerstate1 = samplerstate(rsampler1,smhnprops1,Float64,Float64)
smhsamplerstate2 = samplerstate(rsampler1,smhnprops1,Float64,Float64)
#carry out all the steps of one iteration
srand(rinitseed1) ; initialize!(smhrunner1,rmodel1,smhsamplerstate1,smhchain1,false)
@test smhsamplerstate1.from.values == smhsamplerstate1.proposals.values
propose!(smhsamplerstate1)
@test smhsamplerstate1.from.values != smhsamplerstate1.proposals.values
geometry!(rmodel1,proposals(smhsamplerstate1))
transitionprobability!(smhindicator1,acceptance!(smhsamplerstate1))
sampleindicator!(smhindicator1)
store!(smhrunner1,smhsamplerstate1,smhindicator1,smhchain1)
#perform an iteration
srand(rinitseed1) ; initialize!(smhrunner1,rmodel1,smhsamplerstate2,smhchain2,false)
iterate!(smhrunner2,rmodel1,smhsamplerstate2,smhindicator2,smhchain2,true)
#compare the values from both steps
GMHRunnersTest.testfrom(smhsamplerstate2,smhsamplerstate1.from.values,smhsamplerstate1.from.loglikelihood)
GMHRunnersTest.testproposals(smhsamplerstate2,smhsamplerstate1.proposals.values,smhsamplerstate1.proposals.loglikelihood)

#test the tuning function
println("=========================")
println("Testing of tune! function")
println("=========================")
accepted!(smhtunerstate1,smhindicator1)
accepted!(smhtunerstate2,smhindicator1)
tvals1 = tune(rtuner1,smhtunerstate1)
tvals2 = tune(rtuner2,smhtunerstate2)
cov1 = cov(smhsamplerstate1.density.distribution)

@test index(smhtunerstate1) == 1
tune!(smhrunner1,smhsamplerstate1,rtuner1,smhtunerstate1)
@test index(smhtunerstate1) == 2
@test cov1 == cov(smhsamplerstate1.density.distribution)

@test index(smhtunerstate2) == 1
tune!(smhrunner1,smhsamplerstate1,rtuner2,smhtunerstate2)
@test index(smhtunerstate2) == 2
@test_approx_eq cov1*tvals2*tvals2 cov(smhsamplerstate1.density.distribution)

#run the whole mcmc
println("========================")
println("Testing of run! function")
println("========================")
srand(rinitseed1) ; smhc1 = run!(smhrunner1,rmodel1,rsampler1,rtuner1)
srand(rinitseed1) ; smhc2 = run!(smhrunner2,rmodel1,rsampler1,rtuner2)
@test numsamples(smhc1) == 4
@test numsamples(smhc2) == 7
@test samples(smhc2)[:,1] == rpriorvals1
@test samples(smhc1) == samples(smhc2)[:,4:end]
@test loglikelihood(smhc1) == loglikelihood(smhc2)[4:end]
@test logprior(smhc1,rmodel1) == logprior(smhc2,rmodel1)[4:end]
@test logposterior(smhc1,rmodel1) == logposterior(smhc2,rmodel1)[4:end]

println("==============================")
println("Testing of SMHRunner completed")
println("==============================")


text:
"""
    Binomial(n,p)

A *Binomial distribution* characterizes the number of successes in a sequence of independent trials. It has two parameters: `n`, the number of trials, and `p`, the probability of success in an individual trial, with the distribution:

```math
P(X = k) = {n \\choose k}p^k(1-p)^{n-k},  \\quad \\text{ for } k = 0,1,2, \\ldots, n.
```

```julia
Binomial()      # Binomial distribution with n = 1 and p = 0.5
Binomial(n)     # Binomial distribution for n trials with success rate p = 0.5
Binomial(n, p)  # Binomial distribution for n trials with success rate p

params(d)       # Get the parameters, i.e. (n, p)
ntrials(d)      # Get the number of trials, i.e. n
succprob(d)     # Get the success rate, i.e. p
failprob(d)     # Get the failure rate, i.e. 1 - p
```

External links:

* [Binomial distribution on Wikipedia](http://en.wikipedia.org/wiki/Binomial_distribution)
"""
struct Binomial{T<:Real} <: DiscreteUnivariateDistribution
    n::Int
    p::T

    Binomial{T}(n, p) where {T <: Real} = new{T}(n, p)
end

function Binomial(n::Integer, p::T) where {T <: Real}
    @check_args(Binomial, n >= zero(n))
    @check_args(Binomial, zero(p) <= p <= one(p))
    Binomial{T}(n, p)
end

Binomial(n::Integer, p::T, ::NoArgCheck) where {T<:Real} = Binomial{T}(n, p)
Binomial(n::Integer, p::Integer) = Binomial(n, float(p))
Binomial(n::Integer) = Binomial(n, 0.5)
Binomial() = Binomial(1, 0.5, NoArgCheck())

@distr_support Binomial 0 d.n

#### Conversions

function convert(::Type{Binomial{T}}, n::Int, p::Real) where T<:Real
    return Binomial(n, T(p))
end
function convert(::Type{Binomial{T}}, d::Binomial{S}) where {T <: Real, S <: Real}
    return Binomial(d.n, T(d.p), NoArgCheck())
end

#### Parameters

ntrials(d::Binomial) = d.n
succprob(d::Binomial) = d.p
failprob(d::Binomial{T}) where {T} = one(T) - d.p

params(d::Binomial) = (d.n, d.p)
@inline partype(::Binomial{T}) where {T<:Real} = T

#### Properties

mean(d::Binomial) = ntrials(d) * succprob(d)
var(d::Binomial) = ntrials(d) * succprob(d) * failprob(d)
function mode(d::Binomial{T}) where T<:Real
    (n, p) = params(d)
    n > 0 ? floor(Int, (n + 1) * d.p) : zero(T)
end
modes(d::Binomial) = Int[mode(d)]

median(d::Binomial) = round(Int,mean(d))

function skewness(d::Binomial)
    n, p1 = params(d)
    p0 = 1 - p1
    (p0 - p1) / sqrt(n * p0 * p1)
end

function kurtosis(d::Binomial)
    n, p = params(d)
    u = p * (1 - p)
    (1 - 6u) / (n * u)
end

function entropy(d::Binomial; approx::Bool=false)
    n, p1 = params(d)
    (p1 == 0 || p1 == 1 || n == 0) && return zero(p1)
    p0 = 1 - p1
    if approx
        return (log(twoπ * n * p0 * p1) + 1) / 2
    else
        lg = log(p1 / p0)
        lp = n * log(p0)
        s = exp(lp) * lp
        for k = 1:n
           lp += log((n - k + 1) / k) + lg
           s += exp(lp) * lp
        end
        return -s
    end
end

#### Evaluation & Sampling

@_delegate_statsfuns Binomial binom n p

function rand(rng::AbstractRNG, d::Binomial)
    p, n = d.p, d.n
    if p <= 0.5
        r = p
    else
        r = 1.0-p
    end
    if r*n <= 10.0
        y = rand(rng, BinomialGeomSampler(n,r))
    else
        y = rand(rng, BinomialTPESampler(n,r))
    end
    p <= 0.5 ? y : n-y
end

struct RecursiveBinomProbEvaluator{T<:Real} <: RecursiveProbabilityEvaluator
    n::Int
    coef::T   # p / (1 - p)
end

RecursiveBinomProbEvaluator(d::Binomial) = RecursiveBinomProbEvaluator(d.n, d.p / (1 - d.p))
nextpdf(s::RecursiveBinomProbEvaluator, pv::T, x::Integer) where T <: Real = ((s.n - x + 1) / x) * s.coef * pv

function Base.broadcast!(::typeof(pdf), r::AbstractArray, d::Binomial, X::UnitRange)
    axes(r) == axes(X) || throw(DimensionMismatch("Inconsistent array dimensions."))

    vl,vr, vfirst, vlast = _pdf_fill_outside!(r, d, X)
    if succprob(d) <= 1/2
        # fill normal
        rpe = RecursiveBinomProbEvaluator(d::Binomial)

        # fill central part: with non-zero pdf
        if vl <= vr
            fm1 = vfirst - 1
            r[vl - fm1] = pv = pdf(d, vl)
            for v = (vl + 1):vr
                r[v - fm1] = pv = nextpdf(rpe, pv, v)
            end
        end
    else
        # fill reversed to avoid 1/0 for d.p==1.
        rpe = RecursiveBinomProbEvaluator(d.n, (1 - d.p) / d.p)

        # fill central part: with non-zero pdf
        if vl <= vr
            fm1 = vfirst - 1
            r[vr - fm1] = pv = pdf(d, vr)
            for v = (vr-1):-1:vl
                r[v - fm1] = pv = nextpdf(rpe, pv, d.n-v)
            end
        end
    end
    return r
end
function Base.broadcast(::typeof(pdf), d::Binomial, X::UnitRange)
    r = similar(Array{promote_type(partype(d), eltype(X))}, axes(X))
    r .= pdf.(Ref(d),X)
end

function mgf(d::Binomial, t::Real)
    n, p = params(d)
    (one(p) - p + p * exp(t)) ^ n
end

function cf(d::Binomial, t::Real)
    n, p = params(d)
    (one(p) - p + p * cis(t)) ^ n
end

#### Fit model

struct BinomialStats <: SufficientStats
    ns::Float64   # the total number of successes
    ne::Float64   # the number of experiments
    n::Int        # the number of trials in each experiment

    BinomialStats(ns::Real, ne::Real, n::Integer) = new(ns, ne, n)
end

function suffstats(::Type{<:Binomial}, n::Integer, x::AbstractArray{T}) where T<:Integer
    ns = zero(T)
    for i = 1:length(x)
        @inbounds xi = x[i]
        0 <= xi <= n || throw(DomainError())
        ns += xi
    end
    BinomialStats(ns, length(x), n)
end

function suffstats(::Type{<:Binomial}, n::Integer, x::AbstractArray{T}, w::AbstractArray{Float64}) where T<:Integer
    ns = 0.
    ne = 0.
    for i = 1:length(x)
        @inbounds xi = x[i]
        @inbounds wi = w[i]
        0 <= xi <= n || throw(DomainError())
        ns += xi * wi
        ne += wi
    end
    BinomialStats(ns, ne, n)
end

const BinomData = Tuple{Int, AbstractArray}

suffstats(::Type{<:Binomial}, data::BinomData) = suffstats(Binomial, data...)
suffstats(::Type{<:Binomial}, data::BinomData, w::AbstractArray{Float64}) = suffstats(Binomial, data..., w)

fit_mle(::Type{<:Binomial}, ss::BinomialStats) = Binomial(ss.n, ss.ns / (ss.ne * ss.n))

fit_mle(::Type{<:Binomial}, n::Integer, x::AbstractArray{T}) where {T<:Integer} = fit_mle(Binomial, suffstats(Binomial, n, x))
fit_mle(::Type{<:Binomial}, n::Integer, x::AbstractArray{T}, w::AbstractArray{Float64}) where {T<:Integer} = fit_mle(Binomial, suffstats(Binomial, n, x, w))
fit_mle(::Type{<:Binomial}, data::BinomData) = fit_mle(Binomial, suffstats(Binomial, data))
fit_mle(::Type{<:Binomial}, data::BinomData, w::AbstractArray{Float64}) = fit_mle(Binomial, suffstats(Binomial, data, w))

fit(::Type{<:Binomial}, data::BinomData) = fit_mle(Binomial, data)
fit(::Type{<:Binomial}, data::BinomData, w::AbstractArray{Float64}) = fit_mle(Binomial, data, w)


text:
#####
##### Weighted Essentially Non-Oscillatory (WENO) fifth-order advection scheme
#####

using OffsetArrays
using Oceananigans.Grids: with_halo, return_metrics
using Oceananigans.Architectures: arch_array, architecture
using KernelAbstractions.Extras.LoopInfo: @unroll
using Adapt
import Base: show

const two_32 = Int32(2)

const ƞ = Int32(2) # WENO exponent
const ε = 1e-6

abstract type SmoothnessStencil end

struct VorticityStencil <:SmoothnessStencil end
struct VelocityStencil <:SmoothnessStencil end

"""
    struct WENO5{FT, XT, YT, ZT, XS, YS, ZS, WF} <: AbstractUpwindBiasedAdvectionScheme{2}

Weighted Essentially Non-Oscillatory (WENO) fifth-order advection scheme.

$(TYPEDFIELDS)
"""
struct WENO5{FT, XT, YT, ZT, XS, YS, ZS, VI, WF} <: AbstractUpwindBiasedAdvectionScheme{2}

    "coefficient for ENO reconstruction on x-faces"
    coeff_xᶠᵃᵃ::XT
    "coefficient for ENO reconstruction on x-centers"
    coeff_xᶜᵃᵃ::XT
    "coefficient for ENO reconstruction on y-faces"
    coeff_yᵃᶠᵃ::YT
    "coefficient for ENO reconstruction on y-centers"
    coeff_yᵃᶜᵃ::YT
    "coefficient for ENO reconstruction on z-faces"
    coeff_zᵃᵃᶠ::ZT
    "coefficient for ENO reconstruction on z-centers"
    coeff_zᵃᵃᶜ::ZT

    "coefficient for WENO smoothness indicators on x-faces"
    smooth_xᶠᵃᵃ::XS
    "coefficient for WENO smoothness indicators on x-centers"
    smooth_xᶜᵃᵃ::XS
    "coefficient for WENO smoothness indicators on y-faces"
    smooth_yᵃᶠᵃ::YS
    "coefficient for WENO smoothness indicators on y-centers"
    smooth_yᵃᶜᵃ::YS
    "coefficient for WENO smoothness indicators on z-faces"
    smooth_zᵃᵃᶠ::ZS
    "coefficient for WENO smoothness indicators on z-centers"
    smooth_zᵃᵃᶜ::ZS

    "coefficient for WENO reconstruction, optimal weights"
    C3₀ :: FT
    C3₁ :: FT
    C3₂ :: FT

    function WENO5{VI, WF}(coeff_xᶠᵃᵃ::XT, coeff_xᶜᵃᵃ::XT,
                           coeff_yᵃᶠᵃ::YT, coeff_yᵃᶜᵃ::YT,
                           coeff_zᵃᵃᶠ::ZT, coeff_zᵃᵃᶜ::ZT,
                           smooth_xᶠᵃᵃ::XS, smooth_xᶜᵃᵃ::XS,
                           smooth_yᵃᶠᵃ::YS, smooth_yᵃᶜᵃ::YS,
                           smooth_zᵃᵃᶠ::ZS, smooth_zᵃᵃᶜ::ZS,
                           C3₀::FT, C3₁::FT, C3₂::FT) where {FT, XT, YT, ZT, XS, YS, ZS, VI, WF}

            return new{FT, XT, YT, ZT, XS, YS, ZS, VI, WF}(coeff_xᶠᵃᵃ, coeff_xᶜᵃᵃ, coeff_yᵃᶠᵃ, coeff_yᵃᶜᵃ, coeff_zᵃᵃᶠ, coeff_zᵃᵃᶜ,
                                                           smooth_xᶠᵃᵃ, smooth_xᶜᵃᵃ, smooth_yᵃᶠᵃ, smooth_yᵃᶜᵃ, smooth_zᵃᵃᶠ, smooth_zᵃᵃᶜ,
                                                           C3₀, C3₁, C3₂)
    end
end

"""
    WENO5([FT = Float64;] grid = nothing, stretched_smoothness = false, zweno = false)

Construct a fifth-order weigthed essentially non-oscillatory advection scheme. The constructor allows
construction of WENO schemes on either uniform or stretched grids.

Keyword arguments
=================

  - `grid`: (defaults to `nothing`)
  - `stretched_smoothness`: When `true` it results in computing the coefficients for the smoothness
    indicators β₀, β₁ and β₂ so that they account for the stretched `grid`. (defaults to `false`)
  - `zweno`: When `true` implement a Z-WENO formulation for the WENO weights calculation. (defaults to
    `false`)

!!! warn "No support for WENO5 on curvilinear grids"
    Currently, WENO 5th-order advection schemes don't work for for curvilinear grids.
    Providing `WENO5(::AbstractCurvilinearGrid)` defaults to uniform setting, i.e.
    `WENO5(::AbstractCurvilinearGrid) = WENO5()`.

Not providing any keyword argument, `WENO5()` defaults to the uniform 5th-order coefficients ("uniform
setting) in all directions, using a JS-WENO formulation.

```jldoctest; filter = [r".*┌ Warning.*", r".*└ @ Oceananigans.*"]
julia> using Oceananigans

julia> WENO5()
┌ Warning: defaulting to uniform WENO scheme with Float64 precision, use WENO5(grid = grid) if this was not intended
└ @ Oceananigans.Advection .../src/Advection/weno_fifth_order.jl:90
WENO5 advection scheme with:
    ├── X regular
    ├── Y regular
    └── Z regular
```

`WENO5(grid = grid)` defaults to uniform interpolation coefficient for each of the grid directions that
is uniform (`typeof(Δc) <: Number`) while it precomputes the ENO coefficients for reconstruction for all
grid directions that are stretched. (After testing "on-the-fly" calculation of coefficients for stretched
directions ended up being way too expensive and, therefore, is not supported.)

```jldoctest
julia> using Oceananigans

julia> grid = RectilinearGrid(size = (3, 4, 5), x = (0, 1), y = (0, 1), z = [-10, -9, -7, -4, -1.5, 0]);

julia> WENO5(grid = grid)
WENO5 advection scheme with:
    ├── X regular
    ├── Y regular
    └── Z stretched
```

`WENO5(grid = grid, stretched_smoothness = true)` behaves similarly to `WENO5(grid = grid)` but,
additionally, it also computes the smoothness indicators coefficients, ``β₀``, ``β₁``, and ``β₂``,
taking into account the stretched dimensions.

`WENO5(zweno = true)` implements a Z-WENO formulation for the WENO weights calculation

Comments
========

All methods have the roughly the same execution speed except for `stretched_smoothness = true` that
requires more memory and is less computationally efficient, especially on GPUs. In addition, it has
not been found to be much impactful on the tested cases. As such, most of the times we urge users
to use `WENO5(grid = grid)`, as this increases accuracy on a stretched mesh  but does decreases
memory utilization (and also results in a slight speed-up).

(The above claims were made after some preliminary tests. Thus, we still users to perform some
benchmarks/checks before performing, e.g., a large simulation on a "weirdly" stretched grid.)

On the other hand, a Z-WENO formulation is *most of the times* beneficial (also in case of a uniform
mesh) with roughly the same performances (just a slight slowdown). The same can be said for the
stretched `WENO5(grid = grid)` formulation in case of stretched grids.

References
==========

Shu, Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic
    Conservation Laws, 1997, NASA/CR-97-206253, ICASE Report No. 97-65

Castro et al, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation
    laws, 2011, Journal of Computational Physics, 230(5), 1766-1792
"""

WENO5(grid, FT::DataType=Float64; kwargs...) = WENO5(FT; grid = grid, kwargs...)

function WENO5(FT::DataType = Float64;
               grid = nothing,
               coeffs = nothing,
               stretched_smoothness = false,
               zweno = true,
               vector_invariant = nothing)

    if !(grid isa Nothing)
        FT = eltype(grid)
    end

    weno_coefficients = compute_stretched_weno_coefficients(grid, stretched_smoothness, FT)

    if coeffs isa Nothing
        C3₀, C3₁, C3₂ = FT.((3/10, 3/5, 1/10))
    else
        C3₀, C3₁, C3₂ = FT.(coeffs)
    end

    VI = typeof(vector_invariant)

    return WENO5{VI, zweno}(weno_coefficients..., C3₀, C3₁, C3₂)
end

function compute_stretched_weno_coefficients(grid, stretched_smoothness, FT)

    rect_metrics = (:xᶠᵃᵃ, :xᶜᵃᵃ, :yᵃᶠᵃ, :yᵃᶜᵃ, :zᵃᵃᶠ, :zᵃᵃᶜ)

    if grid isa Nothing
        @warn "defaulting to uniform WENO scheme with $(FT) precision, use WENO5(grid = grid) if this was not intended"
        for metric in rect_metrics
            @eval $(Symbol(:coeff_ , metric)) = nothing
            @eval $(Symbol(:smooth_, metric)) = nothing
        end
    else
        !(grid isa RectilinearGrid) && (@warn "WENO on a curvilinear stretched coordinate is not validated, use at your own risk!!")

        metrics      = return_metrics(grid)
        dirsize      = (:Nx, :Nx, :Ny, :Ny, :Nz, :Nz)

        arch     = architecture(grid)
        new_grid = with_halo((4, 4, 4), grid)

        for (dir, metric, rect_metric) in zip(dirsize, metrics, rect_metrics)
            @eval $(Symbol(:coeff_ , rect_metric)) = calc_interpolating_coefficients($FT, $new_grid.$metric, $arch, $new_grid.$dir)
            @eval $(Symbol(:smooth_, rect_metric)) = calc_smoothness_coefficients($FT, $Val($stretched_smoothness), $new_grid.$metric, $arch, $new_grid.$dir)
        end
    end

    return (coeff_xᶠᵃᵃ , coeff_xᶜᵃᵃ , coeff_yᵃᶠᵃ , coeff_yᵃᶜᵃ , coeff_zᵃᵃᶠ , coeff_zᵃᵃᶜ ,
            smooth_xᶠᵃᵃ, smooth_xᶜᵃᵃ, smooth_yᵃᶠᵃ, smooth_yᵃᶜᵃ, smooth_zᵃᵃᶠ, smooth_zᵃᵃᶜ)
end

# Flavours of WENO
const ZWENO                   = WENO5{<:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, <:Any, true}

const WENOVectorInvariantVel{FT, XT, YT, ZT, XS, YS, ZS, VI, WF}  =
      WENO5{FT, XT, YT, ZT, XS, YS, ZS, VI, WF} where {FT, XT, YT, ZT, XS, YS, ZS, VI<:VelocityStencil, WF}
const WENOVectorInvariantVort{FT, XT, YT, ZT, XS, YS, ZS, VI, WF} =
      WENO5{FT, XT, YT, ZT, XS, YS, ZS, VI, WF} where {FT, XT, YT, ZT, XS, YS, ZS, VI<:VorticityStencil, WF}

const WENOVectorInvariant = WENO5{FT, XT, YT, ZT, XS, YS, ZS, VI, WF} where {FT, XT, YT, ZT, XS, YS, ZS, VI<:SmoothnessStencil, WF}

function Base.show(io::IO, a::WENO5{FT, RX, RY, RZ}) where {FT, RX, RY, RZ}
    print(io, "WENO5 advection scheme with: \n",
              "    ├── X $(RX == Nothing ? "regular" : "stretched") \n",
              "    ├── Y $(RY == Nothing ? "regular" : "stretched") \n",
              "    └── Z $(RZ == Nothing ? "regular" : "stretched")" )
end

Adapt.adapt_structure(to, scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS, VI, WF}) where {FT, XT, YT, ZT, XS, YS, ZS, VI, WF} =
     WENO5{VI, WF}(Adapt.adapt(to, scheme.coeff_xᶠᵃᵃ), Adapt.adapt(to, scheme.coeff_xᶜᵃᵃ),
                   Adapt.adapt(to, scheme.coeff_yᵃᶠᵃ), Adapt.adapt(to, scheme.coeff_yᵃᶜᵃ),
                   Adapt.adapt(to, scheme.coeff_zᵃᵃᶠ), Adapt.adapt(to, scheme.coeff_zᵃᵃᶜ),
                   Adapt.adapt(to, scheme.smooth_xᶠᵃᵃ), Adapt.adapt(to, scheme.smooth_xᶜᵃᵃ),
                   Adapt.adapt(to, scheme.smooth_yᵃᶠᵃ), Adapt.adapt(to, scheme.smooth_yᵃᶜᵃ),
                   Adapt.adapt(to, scheme.smooth_zᵃᵃᶠ), Adapt.adapt(to, scheme.smooth_zᵃᵃᶜ),
                   scheme.C3₀, scheme.C3₁, scheme.C3₂)

@inline boundary_buffer(::WENO5) = 2

@inline symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, centered_fourth_order, c)
@inline symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, centered_fourth_order, c)
@inline symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, centered_fourth_order, c)

@inline symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, ::WENO5, u) = symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, centered_fourth_order, u)
@inline symmetric_interpolate_yᵃᶜᵃ(i, j, k, grid, ::WENO5, v) = symmetric_interpolate_yᵃᶜᵃ(i, j, k, grid, centered_fourth_order, v)
@inline symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, ::WENO5, w) = symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, centered_fourth_order, w)

# Unroll the functions to pass the coordinates in case of a stretched grid
@inline left_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme::WENO5, ψ, args...)  = weno_left_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, ψ, i, Face, args...)
@inline left_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme::WENO5, ψ, args...)  = weno_left_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, ψ, j, Face, args...)
@inline left_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme::WENO5, ψ, args...)  = weno_left_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, ψ, k, Face, args...)

@inline right_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme::WENO5, ψ, args...) = weno_right_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, ψ, i, Face, args...)
@inline right_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme::WENO5, ψ, args...) = weno_right_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, ψ, j, Face, args...)
@inline right_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme::WENO5, ψ, args...) = weno_right_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, ψ, k, Face, args...)

@inline left_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme::WENO5, ψ, args...)  = weno_left_biased_interpolate_xᶠᵃᵃ(i+1, j, k, grid, scheme, ψ, i, Center, args...)
@inline left_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme::WENO5, ψ, args...)  = weno_left_biased_interpolate_yᵃᶠᵃ(i, j+1, k, grid, scheme, ψ, j, Center, args...)
@inline left_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::WENO5, ψ, args...)  = weno_left_biased_interpolate_zᵃᵃᶠ(i, j, k+1, grid, scheme, ψ, k, Center, args...)

@inline right_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme::WENO5, ψ, args...) = weno_right_biased_interpolate_xᶠᵃᵃ(i+1, j, k, grid, scheme, ψ, i, Center, args...)
@inline right_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme::WENO5, ψ, args...) = weno_right_biased_interpolate_yᵃᶠᵃ(i, j+1, k, grid, scheme, ψ, j, Center, args...)
@inline right_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::WENO5, ψ, args...) = weno_right_biased_interpolate_zᵃᵃᶠ(i, j, k+1, grid, scheme, ψ, k, Center, args...)

# Stencil to calculate the stretched WENO weights and smoothness indicators
@inline left_stencil_x(i, j, k, ψ, args...) = @inbounds ( (ψ[i-3, j, k], ψ[i-2, j, k], ψ[i-1, j, k]), (ψ[i-2, j, k], ψ[i-1, j, k], ψ[i, j, k]), (ψ[i-1, j, k], ψ[i, j, k], ψ[i+1, j, k]) )
@inline left_stencil_y(i, j, k, ψ, args...) = @inbounds ( (ψ[i, j-3, k], ψ[i, j-2, k], ψ[i, j-1, k]), (ψ[i, j-2, k], ψ[i, j-1, k], ψ[i, j, k]), (ψ[i, j-1, k], ψ[i, j, k], ψ[i, j+1, k]) )
@inline left_stencil_z(i, j, k, ψ, args...) = @inbounds ( (ψ[i, j, k-3], ψ[i, j, k-2], ψ[i, j, k-1]), (ψ[i, j, k-2], ψ[i, j, k-1], ψ[i, j, k]), (ψ[i, j, k-1], ψ[i, j, k], ψ[i, j, k+1]) )

@inline right_stencil_x(i, j, k, ψ, args...) = @inbounds ( (ψ[i-2, j, k], ψ[i-1, j, k], ψ[i, j, k]), (ψ[i-1, j, k], ψ[i, j, k], ψ[i+1, j, k]), (ψ[i, j, k], ψ[i+1, j, k], ψ[i+2, j, k]) )
@inline right_stencil_y(i, j, k, ψ, args...) = @inbounds ( (ψ[i, j-2, k], ψ[i, j-1, k], ψ[i, j, k]), (ψ[i, j-1, k], ψ[i, j, k], ψ[i, j+1, k]), (ψ[i, j, k], ψ[i, j+1, k], ψ[i, j+2, k]) )
@inline right_stencil_z(i, j, k, ψ, args...) = @inbounds ( (ψ[i, j, k-2], ψ[i, j, k-1], ψ[i, j, k]), (ψ[i, j, k-1], ψ[i, j, k], ψ[i, j, k+1]), (ψ[i, j, k], ψ[i, j, k+1], ψ[i, j, k+2]) )

@inline left_stencil_x(i, j, k, ψ::Function, args...) = @inbounds ( (ψ(i-3, j, k, args...), ψ(i-2, j, k, args...), ψ(i-1, j, k, args...)), (ψ(i-2, j, k, args...), ψ(i-1, j, k, args...), ψ(i, j, k, args...)), (ψ(i-1, j, k, args...), ψ(i, j, k, args...), ψ(i+1, j, k, args...)) )
@inline left_stencil_y(i, j, k, ψ::Function, args...) = @inbounds ( (ψ(i, j-3, k, args...), ψ(i, j-2, k, args...), ψ(i, j-1, k, args...)), (ψ(i, j-2, k, args...), ψ(i, j-1, k, args...), ψ(i, j, k, args...)), (ψ(i, j-1, k, args...), ψ(i, j, k, args...), ψ(i, j+1, k, args...)) )
@inline left_stencil_z(i, j, k, ψ::Function, args...) = @inbounds ( (ψ(i, j, k-3, args...), ψ(i, j, k-2, args...), ψ(i, j, k-1, args...)), (ψ(i, j, k-2, args...), ψ(i, j, k-1, args...), ψ(i, j, k, args...)), (ψ(i, j, k-1, args...), ψ(i, j, k, args...), ψ(i, j, k+1, args...)) )

@inline right_stencil_x(i, j, k, ψ::Function, args...) = @inbounds ( (ψ(i-2, j, k, args...), ψ(i-1, j, k, args...), ψ(i, j, k, args...)), (ψ(i-1, j, k, args...), ψ(i, j, k, args...), ψ(i+1, j, k, args...)), (ψ(i, j, k, args...), ψ(i+1, j, k, args...), ψ(i+2, j, k, args...)) )
@inline right_stencil_y(i, j, k, ψ::Function, args...) = @inbounds ( (ψ(i, j-2, k, args...), ψ(i, j-1, k, args...), ψ(i, j, k, args...)), (ψ(i, j-1, k, args...), ψ(i, j, k, args...), ψ(i, j+1, k, args...)), (ψ(i, j, k, args...), ψ(i, j+1, k, args...), ψ(i, j+2, k, args...)) )
@inline right_stencil_z(i, j, k, ψ::Function, args...) = @inbounds ( (ψ(i, j, k-2, args...), ψ(i, j, k-1, args...), ψ(i, j, k, args...)), (ψ(i, j, k-1, args...), ψ(i, j, k, args...), ψ(i, j, k+1, args...)), (ψ(i, j, k, args...), ψ(i, j, k+1, args...), ψ(i, j, k+2, args...)) )

# Stencil for vector invariant calculation of smoothness indicators in the horizontal direction

# Parallel to the interpolation direction! (same as left/right stencil)
@inline tangential_left_stencil_u(i, j, k, ::Val{1}, u)  = @inbounds left_stencil_x(i, j, k, ℑyᵃᶠᵃ, u)
@inline tangential_left_stencil_u(i, j, k, ::Val{2}, u)  = @inbounds left_stencil_y(i, j, k, ℑyᵃᶠᵃ, u)
@inline tangential_left_stencil_v(i, j, k, ::Val{1}, v)  = @inbounds left_stencil_x(i, j, k, ℑxᶠᵃᵃ, v)
@inline tangential_left_stencil_v(i, j, k, ::Val{2}, v)  = @inbounds left_stencil_y(i, j, k, ℑxᶠᵃᵃ, v)

@inline tangential_right_stencil_u(i, j, k, ::Val{1}, u)  = @inbounds right_stencil_x(i, j, k, ℑyᵃᶠᵃ, u)
@inline tangential_right_stencil_u(i, j, k, ::Val{2}, u)  = @inbounds right_stencil_y(i, j, k, ℑyᵃᶠᵃ, u)
@inline tangential_right_stencil_v(i, j, k, ::Val{1}, v)  = @inbounds right_stencil_x(i, j, k, ℑxᶠᵃᵃ, v)
@inline tangential_right_stencil_v(i, j, k, ::Val{2}, v)  = @inbounds right_stencil_y(i, j, k, ℑxᶠᵃᵃ, v)

#####
##### biased pₖ for û calculation
#####

@inline left_biased_p₀(scheme, ψ, args...) = @inbounds sum(coeff_left_p₀(scheme, args...) .* ψ)
@inline left_biased_p₁(scheme, ψ, args...) = @inbounds sum(coeff_left_p₁(scheme, args...) .* ψ)
@inline left_biased_p₂(scheme, ψ, args...) = @inbounds sum(coeff_left_p₂(scheme, args...) .* ψ)

@inline right_biased_p₀(scheme, ψ, args...) = @inbounds sum(coeff_right_p₀(scheme, args...) .* ψ)
@inline right_biased_p₁(scheme, ψ, args...) = @inbounds sum(coeff_right_p₁(scheme, args...) .* ψ)
@inline right_biased_p₂(scheme, ψ, args...) = @inbounds sum(coeff_right_p₂(scheme, args...) .* ψ)

#####
##### Jiang & Shu (1996) WENO smoothness indicators. See also Equation 2.63 in Shu (1998)
#####

@inline left_biased_β₀(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * (3ψ[1] - 4ψ[2] +  ψ[3])^two_32
@inline left_biased_β₁(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * ( ψ[1]         -  ψ[3])^two_32
@inline left_biased_β₂(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * ( ψ[1] - 4ψ[2] + 3ψ[3])^two_32

@inline right_biased_β₀(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * ( ψ[1] - 4ψ[2] + 3ψ[3])^two_32
@inline right_biased_β₁(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * ( ψ[1]         -  ψ[3])^two_32
@inline right_biased_β₂(FT, ψ, ::Type{Nothing}, scheme, args...) = @inbounds FT(13/12) * (ψ[1] - 2ψ[2] + ψ[3])^two_32 + FT(1/4) * (3ψ[1] - 4ψ[2] +  ψ[3])^two_32

#####
##### Stretched smoothness indicators gathered from precomputed values.
##### The stretched values for β coefficients are calculated from
##### Shu, NASA/CR-97-206253, ICASE Report No. 97-65
##### by hardcoding that p(x) is a 2nd order polynomial
#####

@inline function biased_left_β(ψ, scheme, r, dir, i, location)
    @inbounds begin
        stencil = retrieve_left_smooth(scheme, r, dir, i, location)
        wᵢᵢ = stencil[1]
        wᵢⱼ = stencil[2]
        result = 0
        @unroll for j = 1:3
            result += ψ[j] * ( wᵢᵢ[j] * ψ[j] + wᵢⱼ[j] * dagger(ψ)[j] )
        end
    end
    return result
end

@inline function biased_right_β(ψ, scheme, r, dir, i, location)
    @inbounds begin
        stencil = retrieve_right_smooth(scheme, r, dir, i, location)
        wᵢᵢ = stencil[1]
        wᵢⱼ = stencil[2]
        result = 0
        @unroll for j = 1:3
            result += ψ[j] * ( wᵢᵢ[j] * ψ[j] + wᵢⱼ[j] * dagger(ψ)[j] )
        end
    end
    return result
end

@inline left_biased_β₀(FT, ψ, T, scheme, args...) = biased_left_β(ψ, scheme, 0, args...)
@inline left_biased_β₁(FT, ψ, T, scheme, args...) = biased_left_β(ψ, scheme, 1, args...)
@inline left_biased_β₂(FT, ψ, T, scheme, args...) = biased_left_β(ψ, scheme, 2, args...)

@inline right_biased_β₀(FT, ψ, T, scheme, args...) = biased_right_β(ψ, scheme, 2, args...)
@inline right_biased_β₁(FT, ψ, T, scheme, args...) = biased_right_β(ψ, scheme, 1, args...)
@inline right_biased_β₂(FT, ψ, T, scheme, args...) = biased_right_β(ψ, scheme, 0, args...)

#####
##### VectorInvariant reconstruction (based on JS or Z) (z-direction Val{3} is different from x- and y-directions)
#####
##### Z-WENO-5 reconstruction (Castro et al: High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws)
#####
##### JS-WENO-5 reconstruction
#####

for (side, coeffs) in zip([:left, :right], ([:C3₀, :C3₁, :C3₂], [:C3₂, :C3₁, :C3₀]))
    biased_weno5_weights = Symbol(side, :_biased_weno5_weights)
    biased_β₀ = Symbol(side, :_biased_β₀)
    biased_β₁ = Symbol(side, :_biased_β₁)
    biased_β₂ = Symbol(side, :_biased_β₂)

    tangential_stencil_u = Symbol(:tangential_, side, :_stencil_u)
    tangential_stencil_v = Symbol(:tangential_, side, :_stencil_v)

    biased_stencil_z = Symbol(side, :_stencil_z)

    @eval begin
        @inline function $biased_weno5_weights(FT, ψₜ, T, scheme, dir, idx, loc, args...)
            ψ₂, ψ₁, ψ₀ = ψₜ
            β₀ = $biased_β₀(FT, ψ₀, T, scheme, dir, idx, loc)
            β₁ = $biased_β₁(FT, ψ₁, T, scheme, dir, idx, loc)
            β₂ = $biased_β₂(FT, ψ₂, T, scheme, dir, idx, loc)

            if scheme isa ZWENO
                τ₅ = abs(β₂ - β₀)
                α₀ = scheme.$(coeffs[1]) * (1 + (τ₅ / (β₀ + FT(ε)))^ƞ)
                α₁ = scheme.$(coeffs[2]) * (1 + (τ₅ / (β₁ + FT(ε)))^ƞ)
                α₂ = scheme.$(coeffs[3]) * (1 + (τ₅ / (β₂ + FT(ε)))^ƞ)
            else
                α₀ = scheme.$(coeffs[1]) / (β₀ + FT(ε))^ƞ
                α₁ = scheme.$(coeffs[2]) / (β₁ + FT(ε))^ƞ
                α₂ = scheme.$(coeffs[3]) / (β₂ + FT(ε))^ƞ
            end

            Σα = α₀ + α₁ + α₂
            w₀ = α₀ / Σα
            w₁ = α₁ / Σα
            w₂ = α₂ / Σα

            return w₀, w₁, w₂
        end

        @inline function $biased_weno5_weights(FT, ijk, T, scheme, dir, idx, loc, ::Type{VelocityStencil}, u, v)
            i, j, k = ijk

            u₂, u₁, u₀ = $tangential_stencil_u(i, j, k, dir, u)
            v₂, v₁, v₀ = $tangential_stencil_v(i, j, k, dir, v)

            βu₀ = $biased_β₀(FT, u₀, T, scheme, Val(2), idx, loc)
            βu₁ = $biased_β₁(FT, u₁, T, scheme, Val(2), idx, loc)
            βu₂ = $biased_β₂(FT, u₂, T, scheme, Val(2), idx, loc)

            βv₀ = $biased_β₀(FT, v₀, T, scheme, Val(1), idx, loc)
            βv₁ = $biased_β₁(FT, v₁, T, scheme, Val(1), idx, loc)
            βv₂ = $biased_β₂(FT, v₂, T, scheme, Val(1), idx, loc)

            β₀ = 0.5*(βu₀ + βv₀)
            β₁ = 0.5*(βu₁ + βv₁)
            β₂ = 0.5*(βu₂ + βv₂)

            if scheme isa ZWENO
                τ₅ = abs(β₂ - β₀)
                α₀ = scheme.$(coeffs[1]) * (1 + (τ₅ / (β₀ + FT(ε)))^ƞ)
                α₁ = scheme.$(coeffs[2]) * (1 + (τ₅ / (β₁ + FT(ε)))^ƞ)
                α₂ = scheme.$(coeffs[3]) * (1 + (τ₅ / (β₂ + FT(ε)))^ƞ)
            else
                α₀ = scheme.$(coeffs[1]) / (β₀ + FT(ε))^ƞ
                α₁ = scheme.$(coeffs[2]) / (β₁ + FT(ε))^ƞ
                α₂ = scheme.$(coeffs[3]) / (β₂ + FT(ε))^ƞ
            end

            Σα = α₀ + α₁ + α₂
            w₀ = α₀ / Σα
            w₁ = α₁ / Σα
            w₂ = α₂ / Σα

            return w₀, w₁, w₂
        end

        @inline function $biased_weno5_weights(FT, ijk, T, scheme, ::Val{3}, idx, loc, ::Type{VelocityStencil}, u)
            i, j, k = ijk

            u₂, u₁, u₀ = $biased_stencil_z(i, j, k, u)

            β₀ = $biased_β₀(FT, u₀, T, scheme, Val(3), idx, loc)
            β₁ = $biased_β₁(FT, u₁, T, scheme, Val(3), idx, loc)
            β₂ = $biased_β₂(FT, u₂, T, scheme, Val(3), idx, loc)

            if scheme isa ZWENO
                τ₅ = abs(β₂ - β₀)
                α₀ = scheme.$(coeffs[1]) * (1 + (τ₅ / (β₀ + FT(ε)))^ƞ)
                α₁ = scheme.$(coeffs[2]) * (1 + (τ₅ / (β₁ + FT(ε)))^ƞ)
                α₂ = scheme.$(coeffs[3]) * (1 + (τ₅ / (β₂ + FT(ε)))^ƞ)
            else
                α₀ = scheme.$(coeffs[1]) / (β₀ + FT(ε))^ƞ
                α₁ = scheme.$(coeffs[2]) / (β₁ + FT(ε))^ƞ
                α₂ = scheme.$(coeffs[3]) / (β₂ + FT(ε))^ƞ
            end

            Σα = α₀ + α₁ + α₂
            w₀ = α₀ / Σα
            w₁ = α₁ / Σα
            w₂ = α₂ / Σα

            return w₀, w₁, w₂
        end
    end
end

#####
##### Biased interpolation functions
#####

pass_stencil(ψ, i, j, k, stencil) = ψ
pass_stencil(ψ, i, j, k, ::Type{VelocityStencil}) = (i, j, k)

for (interp, dir, val, cT, cS) in zip([:xᶠᵃᵃ, :yᵃᶠᵃ, :zᵃᵃᶠ], [:x, :y, :z], [1, 2, 3], [:XT, :YT, :ZT], [:XS, :YS, :ZS])
    for side in (:left, :right)
        interpolate_func = Symbol(:weno_, side, :_biased_interpolate_, interp)
        stencil       = Symbol(side, :_stencil_, dir)
        weno5_weights = Symbol(side, :_biased_weno5_weights)
        biased_p₀ = Symbol(side, :_biased_p₀)
        biased_p₁ = Symbol(side, :_biased_p₁)
        biased_p₂ = Symbol(side, :_biased_p₂)

        @eval begin
            @inline function $interpolate_func(i, j, k, grid,
                                               scheme::WENO5{FT, XT, YT, ZT, XS, YS, ZS},
                                               ψ, idx, loc, args...) where {FT, XT, YT, ZT, XS, YS, ZS}

                ψ₂, ψ₁, ψ₀ = ψₜ = $stencil(i, j, k, ψ, grid, args...)
                w₀, w₁, w₂ = $weno5_weights(FT, pass_stencil(ψₜ, i, j, k, Nothing), $cS, scheme, Val($val), idx, loc, Nothing, args...)
                return w₀ * $biased_p₀(scheme, ψ₀, $cT, Val($val), idx, loc) +
                       w₁ * $biased_p₁(scheme, ψ₁, $cT, Val($val), idx, loc) +
                       w₂ * $biased_p₂(scheme, ψ₂, $cT, Val($val), idx, loc)
            end

            @inline function $interpolate_func(i, j, k, grid,
                                               scheme::WENOVectorInvariant{FT, XT, YT, ZT, XS, YS, ZS},
                                               ψ, idx, loc, VI, args...) where {FT, XT, YT, ZT, XS, YS, ZS}

                ψ₂, ψ₁, ψ₀ = ψₜ = $stencil(i, j, k, ψ, grid, args...)
                w₀, w₁, w₂ = $weno5_weights(FT, pass_stencil(ψₜ, i, j, k, VI), $cS, scheme, Val($val), idx, loc, VI, args...)
                return w₀ * $biased_p₀(scheme, ψ₀, $cT, Val($val), idx, loc) +
                       w₁ * $biased_p₁(scheme, ψ₁, $cT, Val($val), idx, loc) +
                       w₂ * $biased_p₂(scheme, ψ₂, $cT, Val($val), idx, loc)
            end
        end
    end
end

#####
##### Coefficients for stretched (and uniform) ENO schemes (see Shu NASA/CR-97-206253, ICASE Report No. 97-65)
#####

@inline coeff_left_p₀(scheme::WENO5{FT}, ::Type{Nothing}, args...) where FT = (  FT(1/3),    FT(5/6), - FT(1/6))
@inline coeff_left_p₁(scheme::WENO5{FT}, ::Type{Nothing}, args...) where FT = (- FT(1/6),    FT(5/6),   FT(1/3))
@inline coeff_left_p₂(scheme::WENO5{FT}, ::Type{Nothing}, args...) where FT = (  FT(1/3),  - FT(7/6),  FT(11/6))

@inline coeff_right_p₀(scheme, ::Type{Nothing}, args...) = reverse(coeff_left_p₂(scheme, Nothing, args...))
@inline coeff_right_p₁(scheme, ::Type{Nothing}, args...) = reverse(coeff_left_p₁(scheme, Nothing, args...))
@inline coeff_right_p₂(scheme, ::Type{Nothing}, args...) = reverse(coeff_left_p₀(scheme, Nothing, args...))

@inline coeff_left_p₀(scheme, T, dir, i, loc) = retrieve_coeff(scheme, 0, dir, i ,loc)
@inline coeff_left_p₁(scheme, T, dir, i, loc) = retrieve_coeff(scheme, 1, dir, i ,loc)
@inline coeff_left_p₂(scheme, T, dir, i, loc) = retrieve_coeff(scheme, 2, dir, i ,loc)

@inline coeff_right_p₀(scheme, T, dir, i, loc) = retrieve_coeff(scheme, -1, dir, i ,loc)
@inline coeff_right_p₁(scheme, T, dir, i, loc) = retrieve_coeff(scheme,  0, dir, i ,loc)
@inline coeff_right_p₂(scheme, T, dir, i, loc) = retrieve_coeff(scheme,  1, dir, i ,loc)

@inline retrieve_coeff(scheme, r, ::Val{1}, i, ::Type{Face})   = scheme.coeff_xᶠᵃᵃ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{1}, i, ::Type{Center}) = scheme.coeff_xᶜᵃᵃ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{2}, i, ::Type{Face})   = scheme.coeff_yᵃᶠᵃ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{2}, i, ::Type{Center}) = scheme.coeff_yᵃᶜᵃ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{3}, i, ::Type{Face})   = scheme.coeff_zᵃᵃᶠ[r+2][i]
@inline retrieve_coeff(scheme, r, ::Val{3}, i, ::Type{Center}) = scheme.coeff_zᵃᵃᶜ[r+2][i]

@inline retrieve_left_smooth(scheme, r, ::Val{1}, i, ::Type{Face})   = scheme.smooth_xᶠᵃᵃ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{1}, i, ::Type{Center}) = scheme.smooth_xᶜᵃᵃ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{2}, i, ::Type{Face})   = scheme.smooth_yᵃᶠᵃ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{2}, i, ::Type{Center}) = scheme.smooth_yᵃᶜᵃ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{3}, i, ::Type{Face})   = scheme.smooth_zᵃᵃᶠ[r+1][i]
@inline retrieve_left_smooth(scheme, r, ::Val{3}, i, ::Type{Center}) = scheme.smooth_zᵃᵃᶜ[r+1][i]

@inline retrieve_right_smooth(scheme, r, ::Val{1}, i, ::Type{Face})   = scheme.smooth_xᶠᵃᵃ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{1}, i, ::Type{Center}) = scheme.smooth_xᶜᵃᵃ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{2}, i, ::Type{Face})   = scheme.smooth_yᵃᶠᵃ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{2}, i, ::Type{Center}) = scheme.smooth_yᵃᶜᵃ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{3}, i, ::Type{Face})   = scheme.smooth_zᵃᵃᶠ[r+4][i]
@inline retrieve_right_smooth(scheme, r, ::Val{3}, i, ::Type{Center}) = scheme.smooth_zᵃᵃᶜ[r+4][i]

@inline calc_interpolating_coefficients(FT, coord::OffsetArray{<:Any, <:Any, <:AbstractRange}, arch, N) = nothing
@inline calc_interpolating_coefficients(FT, coord::AbstractRange, arch, N)                              = nothing

@inline calc_smoothness_coefficients(FT, ::Val{false}, args...) = nothing
@inline calc_smoothness_coefficients(FT, ::Val{true}, coord::OffsetArray{<:Any, <:Any, <:AbstractRange}, arch, N) = nothing
@inline calc_smoothness_coefficients(FT, ::Val{true}, coord::AbstractRange, arch, N) = nothing

function calc_interpolating_coefficients(FT, coord, arch, N)

    cpu_coord = Array(parent(coord))
    cpu_coord = OffsetArray(cpu_coord, coord.offsets[1])

    s1 = create_interp_coefficients(FT,-1, cpu_coord, arch, N)
    s2 = create_interp_coefficients(FT, 0, cpu_coord, arch, N)
    s3 = create_interp_coefficients(FT, 1, cpu_coord, arch, N)
    s4 = create_interp_coefficients(FT, 2, cpu_coord, arch, N)

    return (s1, s2, s3, s4)
end

function create_interp_coefficients(FT, r, cpu_coord, arch, N)

    stencil = NTuple{3, FT}[]
    @inbounds begin
        for i = 0:N+1
            push!(stencil, interp_weights(r, cpu_coord, i, 0, -))
        end
    end
    return OffsetArray(arch_array(arch, stencil), -1)
end

function calc_smoothness_coefficients(FT, beta, coord, arch, N)

    cpu_coord = Array(parent(coord))
    cpu_coord = OffsetArray(cpu_coord, coord.offsets[1])

    s1 = create_smoothness_coefficients(FT, 0, -, cpu_coord, arch, N)
    s2 = create_smoothness_coefficients(FT, 1, -, cpu_coord, arch, N)
    s3 = create_smoothness_coefficients(FT, 2, -, cpu_coord, arch, N)
    s4 = create_smoothness_coefficients(FT, 0, +, cpu_coord, arch, N)
    s5 = create_smoothness_coefficients(FT, 1, +, cpu_coord, arch, N)
    s6 = create_smoothness_coefficients(FT, 2, +, cpu_coord, arch, N)

    return (s1, s2, s3, s4, s5, s6)
end

function create_smoothness_coefficients(FT, r, op, cpu_coord, arch, N)

    # derivation written on overleaf

    stencil = NTuple{2, NTuple{3, FT}}[]
    @inbounds begin
        for i = 0:N+1

            bias1 = Int(op == +)
            bias2 = bias1 - 1

            Δcᵢ = cpu_coord[i + bias1] - cpu_coord[i + bias2]

            Bᵢ  = prim_interp_weights(r, cpu_coord, i, bias1, op)
            bᵢ  =      interp_weights(r, cpu_coord, i, bias1, op)
            bₓᵢ = der1_interp_weights(r, cpu_coord, i, bias1, op)
            Aᵢ  = prim_interp_weights(r, cpu_coord, i, bias2, op)
            aᵢ  =      interp_weights(r, cpu_coord, i, bias2, op)
            aₓᵢ = der1_interp_weights(r, cpu_coord, i, bias2, op)

            pₓₓ = der2_interp_weights(r, cpu_coord, i, op)
            Pᵢ  =  (Bᵢ .- Aᵢ)

            wᵢᵢ = Δcᵢ  .* (bᵢ .* bₓᵢ .- aᵢ .* aₓᵢ .- pₓₓ .* Pᵢ)  .+ Δcᵢ^4 .* (pₓₓ .* pₓₓ)
            wᵢⱼ = Δcᵢ  .* (star(bᵢ, bₓᵢ)  .- star(aᵢ, aₓᵢ) .- star(pₓₓ, Pᵢ)) .+
                                                 Δcᵢ^4 .* star(pₓₓ, pₓₓ)

            push!(stencil, (wᵢᵢ, wᵢⱼ))
        end
    end

    return OffsetArray(arch_array(arch, stencil), -1)
end

@inline dagger(ψ)    = (ψ[2], ψ[3], ψ[1])
@inline star(ψ₁, ψ₂) = (ψ₁ .* dagger(ψ₂) .+ dagger(ψ₁) .* ψ₂)

# Integral of ENO coefficients for 2nd order polynomial reconstruction at the face
function prim_interp_weights(r, coord, i, bias, op)

    coeff = ()
    for j = 0:2
        c = 0
        @inbounds begin
            for m = j+1:3
                num = 0
                for l = 0:3
                    if l != m
                        prod = 1
                        sum  = 0
                        for q = 0:3
                            if q != m && q != l
                                prod *= coord[op(i, r-q+1)]
                                sum  += coord[op(i, r-q+1)]
                            end
                        end
                        num += coord[i+bias]^3 / 3 - sum * coord[i+bias]^2 / 2 + prod * coord[i+bias]
                    end
                end
                den = 1
                for l = 0:3
                    if l!= m
                        den *= (coord[op(i, r-m+1)] - coord[op(i, r-l+1)])
                    end
                end
                c += num / den
            end
        end
        coeff = (coeff..., c * (coord[op(i, r-j)] - coord[op(i, r-j+1)]))
    end

    return coeff
end

# Second derivative of ENO coefficients for 2nd order polynomial reconstruction at the face
function der2_interp_weights(r, coord, i, op)

    coeff = ()
    for j = 0:2
        c = 0
        @inbounds begin
            for m = j+1:3
                num = 0
                for l = 0:3
                    if l != m
                        num += 2
                    end
                end
                den = 1
                for l = 0:3
                    if l!= m
                        den *= (coord[op(i, r-m+1)] - coord[op(i, r-l+1)])
                    end
                end
                c += num / den
            end
        end
        coeff = (coeff..., c * (coord[op(i, r-j)] - coord[op(i, r-j+1)]))
    end

    return coeff
end

# first derivative of ENO coefficients for 2nd order polynomial reconstruction at the face
function der1_interp_weights(r, coord, i, bias, op)

    coeff = ()
    for j = 0:2
        c = 0
        @inbounds begin
            for m = j+1:3
                num = 0
                for l = 0:3
                    if l != m
                        sum = 0
                        for q = 0:3
                            if q != m && q != l
                                sum += coord[op(i, r-q+1)]
                            end
                        end
                        num += 2 * coord[i+bias] - sum
                    end
                end
                den = 1
                for l = 0:3
                    if l!= m
                        den *= (coord[op(i, r-m+1)] - coord[op(i, r-l+1)])
                    end
                end
                c += num / den
            end
        end
        coeff = (coeff..., c * (coord[op(i, r-j)] - coord[op(i, r-j+1)]))
    end

    return coeff
end

# ENO coefficients for 2nd order polynomial reconstruction at the face
function interp_weights(r, coord, i, bias, op)

    coeff = ()
    for j = 0:2
        c = 0
        @inbounds begin
            for m = j+1:3
                num = 0
                for l = 0:3
                    if l != m
                        prod = 1
                        for q = 0:3
                            if q != m && q != l
                                prod *= (coord[i+bias] - coord[op(i, r-q+1)])
                            end
                        end
                        num += prod
                    end
                end
                den = 1
                for l = 0:3
                    if l!= m
                        den *= (coord[op(i, r-m+1)] - coord[op(i, r-l+1)])
                    end
                end
                c += num / den
            end
        end
        coeff = (coeff..., c * (coord[op(i, r-j)] - coord[op(i, r-j+1)]))
    end

    return coeff
end


text:
#
# Copyright (c) 2021 Frederic Bruder, Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

import Downloads
using ZipFile

"""
    list_models()

Shows which models are available. Prints both their IDs and their names.
"""
function list_models()
    @info "Models available in package \"$(@__MODULE__)\""
    header = "ID: model name"
    println("\t$(header)")
    println("\t$("-"^length(header))")
    for i in 1:length(modelNames)
        println("\t$i: $(modelNames[i])")
    end
end

"""
    get_model_filename(modelName::Union{AbstractString, Integer}, tool::AbstractString, version::AbstractString)

Get the filename of a model. `modelName` can be a name or ID. Use `list_models()` to show your options.
"""
function get_model_filename(modelName::AbstractString, tool::AbstractString, version::AbstractString, fmiversion::AbstractString="2.0")

    # workaround to use reference FMUs from the Modelica-Repository
    if tool == "ModelicaReferenceFMUs"
        return download_reference_FMU(modelName, version, fmiversion)
    end

    p_tool = joinpath(p_model_bin, tool)

    if !isdir(p_tool)
        g = glob("*", p_model_bin)
        toolnames = map(g) do x
            splitpath(x)[end]
        end
        error("\"$(tool)\" does not specify an existing tool! Pick one of these: $(toolnames)")
    end

    p_ver = joinpath(p_tool, version)

    if !isdir(p_ver)
        g = glob("*", p_tool)
        versions = map(g) do x
            splitpath(x)[end]
        end
        error("\"$(version)\" does not specify an existing version. Pick one of these: $(versions)")
    end

    p_fmiver = joinpath(p_ver, fmiversion)

    if !isdir(p_fmiver)
        g = glob("*", p_ver)
        fmiversions = map(g) do x
            splitpath(x)[end]
        end
        error("\"$(fmiversion)\" does not specify an existing version. Pick one of these: $(fmiversions)")
    end

    p_model = joinpath(p_fmiver, modelName * ".fmu")
    if !isfile(p_model)
        println(p_model)
        error("\"$(modelName)\" does not specify an existing model. Do `list_models()` to list your options.")
    end

    return p_model
end

function get_model_filename(modelID::Integer, tool::AbstractString, version::AbstractString)
    if modelID < 1
        error("Pass a modelID >= 1.")
    end

    if modelID > length(modelNames)
        error("There are only $(length(modelNames)) models to choose from. Pass modelID between 1 and $(length(modelNames)). Do `list_models()` to list your options.")
    end

    return get_model_filename(modelNames[modelID], tool, version)
end

# TODO prevent this function from downloading on EVERY call. replace `print` by a more appropriate logging function or even an assertion? what happens if the requested model is not found? any information for the user?
"""
    download_reference_FMU(modelName::AbstractString, version::AbstractString="0.0.14", fmiversion::AbstractString="2.0")

Returns the path to an FMU from https://github.com/modelica/Reference-FMUs. If neccessary, all available FMUs are downloaded and saved to a temporary directory at first.
"""
function download_reference_FMU(modelName::AbstractString, version::AbstractString="0.0.14", fmiversion::AbstractString="2.0")

    zipPath = Downloads.download("https://github.com/modelica/Reference-FMUs/releases/download/v$(version)/Reference-FMUs-$(version).zip")
    dir = dirname(zipPath)
    zarchive = ZipFile.Reader(zipPath)

    path = joinpath(dir, "$(modelName)/")
    pathToFmu = joinpath(path, "$(modelName).fmu")

    for f in zarchive.files
        if f.name == "$(fmiversion)/$(modelName).fmu"
            if !ispath(path)
                mkdir(path)
            end

            numBytes = write(pathToFmu, read(f))
            if numBytes == 0
                print("Not able to read!")
            end
        end
    end
    close(zarchive)

    return pathToFmu
end

"""
    generate_mos_scripts()
Generate .mos scripts in `$(p_mos_scripts)` to automate the tool-dependent creation of FMUs.
"""
function generate_mos_scripts(; verbose=true)
    for (name, func) in mosGenerators.generators
        open(io -> write(io, func()), joinpath(p_mos_scripts, "$(name).mos"), "w")
    end

    verbose && @info "Generated all mos scripts in $(p_mos_scripts).\nYou can now copy a path to one of those scripts depending on your Modelica tool and have it executed there to generate all model FMUs into $(p_model_src).\nWhen that's done, call `collect_fmus`."
end

"""
    collect_fmus([dst])
Extracts all FMUs found in directory $(p_model_src) into directory `dst` if specified. Otherwise, the FMUs are moved into a temporary directory.
"""
function collect_fmus(p_dst::Union{AbstractString, Nothing}=nothing)

    if isnothing(p_dst)
        _p_dst = mktempdir(cleanup=false)
    else
        _p_dst = p_dst
    end

    fmuPaths = glob("*.fmu", FMIZoo.p_model_src)

    @assert length(fmuPaths) > 0 "Could not find any FMUs in $(p_model_src). Did you run `FMIZoo.generate_mos_scripts` and have a fitting script executed by your Modelica tool?"

    for fmup in fmuPaths
        fn = splitpath(fmup)[end]
        mv(fmup, joinpath(_p_dst, fn))
    end

    @info "Moved $(length(fmuPaths)) FMU file(s) to $(_p_dst)."
    isnothing(p_dst) && @info "Move the contents into the desired folder!"

end


text:
import Base.normalize

normalize(x::Var) = forward0(normalize, x)

function forward{T}(::typeof(normalize), x::Matrix{T})
    z = mapreducedim(v -> v*v, +, x, 1)
    @inbounds @simd for i = 1:length(z)
        z[i] = T(1) / sqrt(z[i]+1e-9)
    end
    y = x .* z
    backward!(gy, gx) = isvoid(gx) || ∇normalize!(gy, x, gx, z)
    y, backward!
end

function ∇normalize!{T}(gy::Matrix{T}, x::Matrix{T}, gx::Matrix{T}, z::Matrix{T})
    s = sum(x, 1)
    for j = 1:size(x,2)
        @inbounds @simd for i = 1:size(x,1)
            gx[i,j] += gy[i,j] * (z[j] - z[j]*z[j]*z[j]*s[j]*x[i,j])
        end
    end
end


text:
"""
    CTask

Wrapper of a [`Task`](@ref) for which deep copying of stack allocated objects is enabled.
"""
struct CTask
    task::Task

    function CTask(task::Task)
        ret = new(enable_stack_copying(task))
        task.storage === nothing && (task.storage = IdDict())
        task.storage[:ctask] = ret
        ret
    end
end

CTask(f) = CTask(Task(task_wrapper(f)))

# Iteration interface.
Base.iterate(ctask::CTask, state=nothing) = consume(ctask), nothing
Base.IteratorSize(::Type{CTask}) = Base.SizeUnknown()
Base.IteratorEltype(::Type{CTask}) = Base.EltypeUnknown()

struct CTaskException <: Exception
    task::Task
end

function Base.showerror(io::IO, ex::CTaskException)
    println(io, "CTaskException:")
    ct = ex.task
    bt = @static if VERSION < v"1.6.0-DEV.1145"
        ct.backtrace
    else
        ct.storage[:_libtask_bt]
    end
    showerror(io, ct.exception, bt)
end

# Utility function for self-copying mechanism
_current_task() = ccall((:vanilla_get_current_task, libtask_julia), Ref{Task}, ())

n_copies() = n_copies(_current_task())
function n_copies(t::Task)
  t.storage === nothing && (t.storage = IdDict())
  return get!(t.storage, :n_copies, 0)
end

function enable_stack_copying(t::Task)
    if istaskfailed(t)
        error("only runnable or finished tasks' stack can be copied.")
    end
    # return ccall((:jl_enable_stack_copying, libtask_julia),
    #             Any, (Any,), t)::Task
    copy_stack = internal_getfield(t, :copy_stack, Int32)
    if copy_stack == 0
        internal_setfield(t, :copy_stack, Int32(1))
        internal_setfield(t, :bufsz, Csize_t(0))
        ccall((:jl_reset_task_ctx, libtask_julia),
              Cvoid, (Any, Csize_t, Csize_t, Csize_t),
              t, TASK_OFFSETS[:ctx], TASK_OFFSETS[:sizeof_ctx],
              TASK_OFFSETS[:tls_base_context])
    end
    return t
end

"""
    task_wrapper()

`task_wrapper` is a wordaround for set the result/exception to the
correct task which maybe copied/forked from another one(the original
one). Without this, the result/exception is always sent to the
original task. That is done in `JULIA_PROJECT/src/task.c`, the
function `start_task` and `finish_task`.

This workaround is not the proper way to do the work it does. The
proper way is refreshing the `current_task` (the variable `t`) in
`start_task` after the call to `jl_apply` returns.

"""
function task_wrapper(func)
    f = let func=func
        () -> begin
            try
                ct = _current_task()
                res = func()
                ct.result = res
                ct.storage === nothing && (ct.storage = IdDict())
                ct.storage[:_libtask_state] = :done
                wait()
            catch ex
                ct = _current_task()
                @static if VERSION < v"1.6.0-DEV.1145"
                    ct.exception = ex
                    ct.backtrace = catch_backtrace()
                else
                    ct._isexception = true
                    ct.storage[:_libtask_bt] = catch_backtrace()
                end
                ct.result = ex
                ct.storage === nothing && (ct.storage = IdDict())
                ct.storage[:_libtask_state] = :failed
                wait()
            end
        end
    end
    return f
end

function Base.copy(ctask::CTask)
    task = ctask.task
    if istaskfailed(task)
        error("only runnable or finished tasks can be copied.")
    end

    # memory copy
    # newtask = ccall((:jl_clone_task, libtask_julia), Any, (Any,), task)::Task
    newtask = ccall((:jl_clone_task_opaque, libtask_julia),
                    Any, (Any, Csize_t), task, TASK_OFFSETS[:END])::Task
    internal_setfield(newtask, :exception, nothing)
    internal_setfield(newtask, :backtrace, nothing)
    internal_setfield(newtask, :tls, nothing)
    internal_setfield(newtask, :result, nothing)
    internal_setfield(newtask, :donenotify, nothing)
    internal_setfield(newtask, :excstack, C_NULL)
    internal_setfield(newtask, :ptls, C_NULL)
    internal_setfield(newtask, :gcstack, C_NULL)

    if haskey(TASK_OFFSETS, :stkbuf) && haskey(TASK_OFFSETS, :bufsz)
        old_stkbuf = internal_getfield(task, :stkbuf, Ptr)
        if old_stkbuf != C_NULL
            internal_setfield(task, :bufsz, Csize_t(0))
        else
            internal_setfield(newtask, :stkbuf, C_NULL)
        end
        internal_setfield(newtask, :bufsz, Csize_t(0))
    end
    memset(newtask, 0, :locks)
    # memory copy done

    task.storage[:n_copies] = 1 + n_copies(task)
    newtask.storage = copy(task.storage)

    # copy fields not accessible in task.c
    newtask.code = task.code
    setstate!(newtask, getstate(task))
    newtask.result = task.result

    copy_tarrays(task, newtask)
    return CTask(newtask)
end

function produce(v)
    ct = _current_task()
    ct.storage === nothing && (ct.storage = IdDict())

    consumers = get!(ct.storage, :consumers, nothing)
    local empty, task
    while true
        if consumers isa Task
            task = consumers
            ct.storage[:consumers] = nothing
            empty = true
            break
        elseif consumers isa Condition && !isempty(consumers.waitq)
            task = popfirst!(consumers.waitq)
            empty = isempty(consumers.waitq)
            break
        end

        # Wait until there are more consumers.
        wait()

        # Update consumers.
        consumers = ct.storage[:consumers]
    end

    # Internal check to make sure that it is possible to switch to the consumer.
    @assert !istaskdone(task) && !istaskfailed(task)

    task.queue !== nothing && yield()

    if empty
        # Switch to the consumer.
        schedule(task, v)
        wait()

        ct = _current_task() # When a task is copied, ct should be updated to new task ID.
        while true
            q = ct.storage[:consumers]
            if isa(q,Task)
                return q.result
            elseif isa(q,Condition) && !isempty(q.waitq)
                return q.waitq[1].result
            end

            # Wait until there are more consumers.
            wait()
        end
    else
        schedule(task, v)
        # make sure `task` runs before us. otherwise, the producer might
        # finish before `t` runs again, causing it to see the producer
        # as done, causing done(::Task, _) to miss the value `v`.
        # see issue #7727
        yield()
        return consumers.waitq[1].result
    end
end

function consume(ctask::CTask, values...)
    # Check if the producer is done.
    producer = ctask.task
    istaskdone(producer) && return wait(producer)

    # Obtain the current task and set its result.
    ct = _current_task()
    ct.result = length(values) == 1 ? values[1] : values

    # Obtain the consumers listening for the producer.
    producer.storage === nothing && (producer.storage = IdDict())
    consumers = get!(producer.storage, :consumers, nothing)

    if consumers === nothing || (consumers isa Condition && isempty(consumers.waitq))
        # Set the current task as consumer if none exists.
        producer.storage[:consumers] = ct
    else
        # Otherwise add the current task to the waiting queue.
        if consumers isa Task
            # If there is no queue currently but only a single task, replace it with a
            # waiting queue with this task.
            producer.storage[:consumers] = Condition()
            push!(producer.storage[:consumers].waitq, consumers)
        end
        push!(producer.storage[:consumers].waitq, ct)
    end

    if !istaskdone(producer) && !istaskfailed(producer)
        # Switch to the producer.
        schedule(producer)
        yield()

        # Update the state if possible.
        if producer.storage isa IdDict && haskey(producer.storage, :_libtask_state)
            setstate!(producer, producer.storage[:_libtask_state])
        end

        # If the task failed, throw an exception.
        istaskfailed(producer) && throw(CTaskException(producer))

        # If the task is done return the result.
        istaskdone(producer) && return producer.result
    end

    wait()
end

function getstate(task::Task)
    @static if VERSION < v"1.6.0-DEV.618"
        return task.state
    else
        return task._state
    end
end

function setstate!(task::Task, state)
    @static if VERSION < v"1.6.0-DEV.618"
        task.state = state
    else
        if state === :runnable
            task._state = Base.task_state_runnable
        elseif state === :done
            task._state = Base.task_state_done
        elseif state === :failed
            task._state = Base.task_state_failed
        else
            task._state = state
        end
    end
end