text:
"""
`module  JAC.DecayYield`
    ... a submodel of JAC that contains all methods for computing fluorescence and Auger yields for some level(s).
"""
module DecayYield

    using Printf, JAC, ..Basics, ..Defaults, ..ManyElectron, ..Radial, ..TableStrings

    """
    `struct  DecayYield.Settings  <:  AbstractPropertySettings`  ... defines a type for the details and parameters of computing fluorescence and Auger yields.

        + approach                 ::String         ... Determines the applied cascade approach for the yields:
                                                        {"AverageSCA", "SCA"}.
        + printBefore              ::Bool           ... True if a list of selected levels is to be printed before the
                                                        actual computations start.
        + levelSelection           ::LevelSelection ... Specifies the selected levels, if any.
    """
    struct Settings  <:  AbstractPropertySettings
        approach                   ::String
        printBefore                ::Bool
        levelSelection             ::LevelSelection
    end

    """
    `DecayYield.Settings()`
        ... constructor for an `empty` instance of DecayYield.Settings for the computation of fluorescence and Auger yields.
    """
    function Settings()
        Settings("AverageSCA", false, LevelSelection() )
    end

    # `Base.show(io::IO, settings::DecayYield.Settings)`  ... prepares a proper printout of the variable settings::DecayYield.Settings.
    function Base.show(io::IO, settings::DecayYield.Settings)
        println(io, "approach:                 $(settings.approach)  ")
        println(io, "printBefore:              $(settings.printBefore)  ")
        println(io, "levelSelection:           $(settings.levelSelection)  ")
    end

    """
    `struct  DecayYield.Outcome`
        ... defines a type to keep the outcome of a fluorescence and Auger yield computation as well other results.

        + level                     ::Level              ... Atomic level to which the outcome refers to.
        + NoPhotonLines             ::Int64              ... Number of photon decay lines
        + NoAugerLines              ::Int64              ... Number of Auger decay lines
        + rateR                     ::EmProperty         ... Total fluorescence rate
        + rateA                     ::Float64            ... Total Auger rate
        + omegaR                    ::EmProperty         ... fluorescence yield
        + omegaA                    ::EmProperty         ... Auger yield
    """
    struct Outcome
        level                       ::Level
        NoPhotonLines               ::Int64
        NoAugerLines                ::Int64
        rateR                       ::EmProperty
        rateA                       ::Float64
        omegaR                      ::EmProperty
        omegaA                      ::EmProperty
    end

    """
    `DecayYield.Outcome()`
        ... constructor for an `empty` instance of DecayYield.Outcome for the computation of fluorescence and Auger yields.
    """
    function Outcome()
        Outcome(Level(), 0, 0, EmProperty(0.), 0., EmProperty(0.), EmProperty(0.))
    end

    # `Base.show(io::IO, outcome::DecayYield.Outcome)`  ... prepares a proper printout of the variable outcome::DecayYield.Outcome.
    function Base.show(io::IO, outcome::DecayYield.Outcome)
        println(io, "level:                $(outcome.level)  ")
        println(io, "NoPhotonLines:        $(outcome.NoPhotonLines)  ")
        println(io, "NoAugerLines:         $(outcome.NoAugerLines)  ")
        println(io, "NoPhotonLines:        $(outcome.NoPhotonLines)  ")
        println(io, "NoPhotonLines:        $(outcome.NoPhotonLines)  ")
    end

    """
    `DecayYield.computeOutcomes(configs::Array{Configuration,1},asfSettings::AsfSettings,
                                multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::DecayYield.Settings; output=true)`
        ... to compute (as selected) the alpha-variation parameters for the levels of the given multiplet and as specified by
            the given settings. The results are printed in neat tables to screen but nothing is returned otherwise.
    """
    function computeOutcomes(initialConfigs::Array{Configuration,1}, asfSettings::AsfSettings,
                             multiplet::Multiplet, nm::Nuclear.Model, grid::Radial.Grid, settings::DecayYield.Settings; output=true)
        println("")
        printstyled("DecayYield.computeOutcomes(): The computation of the fluorescence & Auger yields starts now ... \n", color=:light_green)
        printstyled("----------------------------------------------------------------------------------------------- \n", color=:light_green)
        #
        outcomes = DecayYield.determineOutcomes(multiplet, settings)
        # Display all selected levels before the computations start
        if  settings.printBefore    DecayYield.displayOutcomes(outcomes)    end
        #
        # Overwrite the default settings for the generation of the continuum orbitals, if needed
        ## setDefaults("method: continuum, asymptotic Coulomb")
        ## setDefaults("method: normalization, pure sine")
        #
        println("\nPerform a cascade computation for all radiative decay channels of the levels from the initial configurations:")
        if      settings.approach == "AverageSCA"    cApproach = Cascade.AverageSCA()
        elseif  settings.approach == "SCA"           cApproach = Cascade.SCA()
        else    error("Improper (cascade) approach for yield computations; approach = $(settings.approach)")
        end
        wa = Cascade.Computation(Cascade.Computation(), name="photon lines", nuclearModel=nm, grid=grid, asfSettings=asfSettings,
                                 scheme=Cascade.StepwiseDecayScheme([Radiative()], 0, Dict{Int64,Float64}(), 0, Shell[], Shell[]),
                                 approach=cApproach, initialConfigs=initialConfigs)
        wb = perform(wa, output=true, outputToFile=false);   linesR = wb["decay line data:"].linesR
        println("\nPerform a cascade computation for all (single-electron) Auger decay channels of the levels from the initial configurations:")
        wa = Cascade.Computation(Cascade.Computation(), name="Auger lines", nuclearModel=nm, grid=grid, asfSettings=asfSettings,
                                 scheme=Cascade.StepwiseDecayScheme([Auger()], 1, Dict{Int64,Float64}(), 0, Shell[], Shell[]),
                                 approach=cApproach, initialConfigs=initialConfigs)
        wb = perform(wa, output=true, outputToFile=false);   linesA = wb["decay line data:"].linesA
        #
        # Calculate all amplitudes and requested properties
        newOutcomes = DecayYield.Outcome[]
        for  outcome in outcomes
            newOutcome = Cascade.computeDecayYieldOutcome(outcome, linesR, linesA, settings)
            push!( newOutcomes, newOutcome)
        end
        # Print all results to screen
        DecayYield.displayResults(stdout, newOutcomes, settings)
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary    DecayYield.displayResults(iostream, newOutcomes, settings)   end
        #
        if    output    return( newOutcomes )
        else            return( nothing )
        end
    end

    """
    `DecayYield.determineOutcomes(multiplet::Multiplet, settings::DecayYield.Settings)`
        ... to determine a list of Outcomes's for the computation of the alpha-variation parameters for the given
            multiplet. It takes into account the particular selections and settings. An Array{DecayYield.Outcome,1} is
            returned. Apart from the level specification, all physical properties are set to zero during the
            initialization process.
    """
    function  determineOutcomes(multiplet::Multiplet, settings::DecayYield.Settings)
        outcomes = DecayYield.Outcome[]
        for  level  in  multiplet.levels
            if  Basics.selectLevel(level, settings.levelSelection)
                push!( outcomes, DecayYield.Outcome(level, 0, 0, EmProperty(0.), 0., EmProperty(0.), EmProperty(0.)) )
            end
        end
        return( outcomes )
    end

    """
    `DecayYield.displayOutcomes(outcomes::Array{DecayYield.Outcome,1})`
        ... to display a list of levels that have been selected for the computations. A small neat table of all
            selected levels and their energies is printed but nothing is returned otherwise.
    """
    function  displayOutcomes(outcomes::Array{DecayYield.Outcome,1})
        nx = 43
        println(" ")
        println("  Selected DecayYield levels:")
        println(" ")
        println("  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=4);
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=4)
        println(sa);    println(sb);    println("  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", outcome.level.energy)) * "    "
            println( sa )
        end
        println("  ", TableStrings.hLine(nx))
        #
        return( nothing )
    end

    """
    `DecayYield.displayResults(stream::IO, outcomes::Array{DecayYield.Outcome,1}, settings::DecayYield.Settings)`
        ... to display the energies, M_ms and F-parameters, etc. for the selected levels. A neat table is printed but nothing
            is returned otherwise.
    """
    function  displayResults(stream::IO, outcomes::Array{DecayYield.Outcome,1}, settings::DecayYield.Settings)
        nx = 149
        println(stream, " ")
        println(stream, "  Fluorescence and Auger decay yields:")
        println(stream, " ")
        println(stream, "    + No_R, rate_R, omega_R  ... number of fluorescence lines, total fluorescence rate and yield. ")
        println(stream, "    + No_A, rate_A, omega_A  ... number of Auger lines, total Auger rate and Auger yield. ")
        if    settings.approach in ["AverageSCA", "SCA"]
        println(stream, "    + Approach:  $(settings.approach)  ... all fluorescence rates and yields only in Babushkin gauge. ")
        end
        #
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  ";   sb = "  "
        sa = sa * TableStrings.center(10, "Level"; na=2);                             sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(10, "J^P";   na=4);                             sb = sb * TableStrings.hBlank(14)
        sa = sa * TableStrings.center(14, "Energy"; na=3)
        sb = sb * TableStrings.center(14, TableStrings.inUnits("energy"); na=3)
        sa = sa * TableStrings.center( 9, "No_R"    ; na=0);                          sb = sb * TableStrings.hBlank(9)
        sa = sa * TableStrings.center(24, "Cou -- rate_R  -- Bab"; na=1);             sb = sb * TableStrings.hBlank(25)
        sa = sa * TableStrings.center(24, "Cou -- omega_R -- Bab"; na=3);             sb = sb * TableStrings.hBlank(27)
        sa = sa * TableStrings.center( 6, "No_A"    ; na=3);                          sb = sb * TableStrings.hBlank(9)
        sa = sa * TableStrings.center(10, "rate_A"; na=2);                            sb = sb * TableStrings.hBlank(12)
        sa = sa * TableStrings.center(24, "Cou -- omega_A -- Bab"; na=3);             sb = sb * TableStrings.hBlank(27)
        println(stream, sa);    println(stream, sb);    println(stream, "  ", TableStrings.hLine(nx))
        #
        for  outcome in outcomes
            sa  = "  ";    sym = LevelSymmetry( outcome.level.J, outcome.level.parity)
            sa = sa * TableStrings.center(10, TableStrings.level(outcome.level.index); na=2)
            sa = sa * TableStrings.center(10, string(sym); na=4)
            energy = outcome.level.energy
            sa = sa * @sprintf("%.8e", Defaults.convertUnits("energy: from atomic", energy))    * "     "
            sa = sa * string(outcome.NoPhotonLines, pad=4)                                      * "    "
            sa = sa * @sprintf("%.4e", outcome.rateR.Coulomb)                                   * " "
            sa = sa * @sprintf("%.4e", outcome.rateR.Babushkin)                                 * "    "
            sa = sa * @sprintf("%.4e", outcome.omegaR.Coulomb)                                  * " "
            sa = sa * @sprintf("%.4e", outcome.omegaR.Babushkin)                                * "     "
            sa = sa * string(outcome.NoAugerLines, pad=4)                                       * "    "
            sa = sa * @sprintf("%.4e", outcome.rateA)                                           * "    "
            sa = sa * @sprintf("%.4e", outcome.omegaA.Coulomb)                                  * " "
            sa = sa * @sprintf("%.4e", outcome.omegaA.Babushkin)                                * "  "
            println(stream, sa )
        end
        println(stream, "  ", TableStrings.hLine(nx), "\n\n")
        #
        return( nothing )
    end

end # module


text:
@doc raw"""
    gaussradau(n::Integer) -> Tuple{Vector{Float64},Vector{Float64}}

Return nodes and weights of [Gauss-Radau quadrature](https://mathworld.wolfram.com/RadauQuadrature.html).

```math
\int_{-1}^{1} f(x) dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gaussradau(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 2/5
true
```

Note that the first node is fixed at -1.

```jldoctest
julia> x, w = gaussradau(3);

julia> x[1]
-1.0
```
"""
function gaussradau(n::Integer, T::Type=Float64)
    a = b = zero(T)
    # RADAUPTS   Gauss-Legendre-Radau Quadrature Nodes and Weights
    if n == 1
        return T[-1], T[2]
    elseif n == 2
        return [-1, one(T)/3], [one(T)/2, convert(T,3)/2]
    else
        # Compute via GaussJacobi:
        x, w = gaussjacobi(n - 1, a, b+1)
        @inbounds for i in 1:length(w)
            w[i] = w[i] / (1 + x[i])
        end
        pushfirst!(x, -1)
        pushfirst!(w, convert(T, 2) / n^2)
        return x, w
    end
end

function gaussradau(n::Integer, α, β)
    if n ≤ 0
        throw(DomainError(n, "Input N must be a positive integer"))
    end
    m = n - 1
    s = α + β
    T = float(eltype(s))
    μ = jacobimoment(α, β)
    n == 0 && return T[], T[]
    n == 1 && return [-one(T)], [μ]
    J = jacobi_jacobimatrix(n, α, β)
    aᴿ = -1 + 2m*convert(T,m+α)/((2m+s)*(2m+s+1))
    J.dv[end] = aᴿ
    x, V = eigen(J)
    w = V[1,:].^2 .* μ
    x[1] = -1 # fix rounding from eigen computation
    return x, w
end

ext: jl


text:
import Ghost: FunctionResolver, __new__

const DRULES = Ref(FunctionResolver{Function}())

function expr2signature(modl::Module, ex::Expr)
    BAD_FORMAT_MSG = "Expected signature as `fn(x::T1, y::T2)`"
    @assert Meta.isexpr(ex, :call) BAD_FORMAT_MSG
    # instead of modl.eval() we could use getfield(modl, :name)
    # but it doesn't work with complex stuff like qualified functions
    # or parametric types
    fn_typ = modl.eval(ex.args[1]) |> typeof
    arg_typs = []
    for x in ex.args[2:end]
        @assert Meta.isexpr(x, :(::)) BAD_FORMAT_MSG
        # Module.eval() is slower, but supports UnionAll
        typ = modl.eval(x.args[end])
        push!(arg_typs, typ)
    end
    return Tuple{fn_typ, arg_typs...}
end

macro drule(ex, df)
    let sig = expr2signature(@__MODULE__, ex)
        quote
            $DRULES[][$sig] = $df
            nothing
        end
    end
end

get_deriv_function(sig) = DRULES[][sig]

is_yota_primitive(sig) = sig in DRULES[]

###############################################################################
#                                    RULES                                    #
###############################################################################

# @diffrule *(u::AbstractArray, v::AbstractArray)    u     dy * transpose(v)

# @diffrule *(u::Real         , v::Real )            v     u * dy
# @diffrule *(u::Real         , v::AbstractArray)    v     u .* dy
# @diffrule *(u::AbstractArray, v::Real )            v     sum(u .* dy)
# @diffrule *(u::AbstractArray, v::AbstractArray)    v     transpose(u) * dy

import Base.Broadcast: broadcasted, materialize

function ∇times(dy, ::typeof(*), x::AbstractArray, y::AbstractArray)
    return NoTangent(), dy * y', x' * dy
end
@drule *(x::AbstractArray, y::AbstractArray) ∇times

function ∇times(dy, ::typeof(*), x::Number, y::Number)
    return NoTangent(), dy * y, dy * x
end
@drule *(x::Number, y::Number) ∇times

function ∇plus(dy, ::typeof(+), x::Number, y::Number)
    return NoTangent(), dy, dy
end
@drule +(x::Number, y::Number) ∇plus

# function ∇broadcasted(dy, ::typeof(Broadcast.broadcasted), f::F, args...) where F
#     df = get_deriv_function(Tuple{F, map(eltype, args)...})
#     # TODO: write a more descriptive error message
#     @assert(df !== nothing,
#         "Broadcasting using Yota rule on a non-primitive function $f")
#     # TODO: use this approach instead:
#     # https://github.com/JuliaDiff/ChainRulesCore.jl/issues/68#issuecomment-834437769
#     darg_tuples = df.(dy, f, args...)
#     darg_single_tuple = tuple(
#         [vcat([t[i] for t in darg_tuples]) for i in 1:length(args) + 1]...
#     )
#     return NoTangent(), darg_single_tuple
# end
# @drule Broadcast.broadcasted(f::Function, args::Vararg) ∇broadcasted

function ∇materialize(dy, ::typeof(Broadcast.materialize), x::Any)
    return NoTangent(), dy
end
@drule Broadcast.materialize(x::Any) ∇materialize

# some special broadcasting rules

# are these rules just incorrect versions of the `∇broadcasted_special` below?
# function ∇broadcasted(dy, ::typeof(Broadcast.broadcasted), f::typeof(+), args...)
#     return NoTangent(), NoTangent(), [dy for a in args]...
# end
# @drule Broadcast.broadcasted(f::typeof(+), args::Vararg) ∇broadcasted

# function ∇broadcasted(dy, ::typeof(Broadcast.broadcasted), f::typeof(*), args...)
#     return NoTangent(), NoTangent(), [dy .* a for a in args]...
# end
# @drule Broadcast.broadcasted(f::typeof(*), args::Vararg) ∇broadcasted

# @diffrule getindex(u::AbstractArray, i)         u    ungetindex(u, dy, i)

function ∇getindex(dy, ::typeof(getindex), x, I...)
    return NoTangent(), ungetindex(x, dy, I...), [ZeroTangent() for i in I]...
end
@drule getindex(x::Any, I::Vararg) ∇getindex

function ∇sum(dy, ::typeof(sum), x::AbstractArray)
    dx = similar(x)
    dx .= dy
    return NoTangent(), dx
end
@drule sum(x::AbstractArray) ∇sum

function ∇mean(dy, ::typeof(Statistics.mean), x::AbstractArray, dims=1:ndims(x))
    dx = similar(x)
    dx .= dy ./ prod(size(x, d) for d in dims)
    return NoTangent(), dx
end
@drule Statistics.mean(x::AbstractArray) ∇mean

## special broadcast

function ∇broadcasted_special(dy, ::typeof(broadcasted), ::typeof(+), x, y)
    return NoTangent(), NoTangent(), unbroadcast(x, dy), unbroadcast(y, dy)
end
@drule broadcasted(::typeof(+), x::AbstractArray, y::AbstractArray) ∇broadcasted_special

function ∇broadcasted_special(dy, ::typeof(broadcasted), ::typeof(*), x, y)
    return NoTangent(), NoTangent(), unbroadcast_prod_x(x, y, dy), unbroadcast_prod_y(x, y, dy)
end
@drule broadcasted(::typeof(*), x::AbstractArray, y::AbstractArray) ∇broadcasted_special

function ∇broadcasted(dy, ::typeof(broadcasted), ::typeof(Base.literal_pow),
    ::typeof(^), x::Any, ::Val{p}) where p
    return NoTangent(), NoTangent(), NoTangent(), (@. p * x ^ (p - 1) * dy), ZeroTangent()
end
@drule broadcasted(::typeof(Base.literal_pow), ::typeof(^), x::Any, ::Val) ∇broadcasted

function ∇broadcasted(dy, ::typeof(broadcasted),
    ::typeof(^), x::Any, p::Real)
    return NoTangent(), NoTangent(), (@. p * x ^ (p - 1) * dy), ZeroTangent()
end
@drule broadcasted(::typeof(^), x::Any, ::Real) ∇broadcasted

## __new__, getfield, getproperty

function ∇getproperty(dy, ::typeof(getproperty), s, f::Symbol)
    T = typeof(s)
    nt = NamedTuple{(f,)}((dy,))
    return NoTangent(), Tangent{T}(; nt...), ZeroTangent()
end
@drule getproperty(s::Any, f::Symbol) ∇getproperty

function ∇getfield(dy, ::typeof(getfield), s::Tuple, f::Int)
    T = typeof(s)
    return NoTangent(), Tangent{T}([i == f ? dy : ZeroTangent() for i=1:length(s)]...), ZeroTangent()
end
@drule getfield(s::Tuple, f::Union{Symbol, Int}) ∇getfield

function ∇__new__(dy, ::typeof(__new__), T, args...)
    if dy isa NoTangent || dy isa ZeroTangent
        fld_derivs = [dy for fld in fieldnames(T)]
    else
        fld_derivs = [getproperty(dy, fld) for fld in fieldnames(T)]
    end
    return NoTangent(), NoTangent(), fld_derivs...
end
@drule __new__(T::Any, args::Vararg) ∇__new__

## iterate

function ∇iterate(dy, ::typeof(iterate), x::AbstractArray)
    return NoTangent(), ungetindex(x, dy, 1)
end
function ∇iterate(dy, ::typeof(iterate), x::AbstractArray, i::Integer)
    return NoTangent(), ungetindex(x, dy, i), ZeroTangent()
end
@drule iterate(x::AbstractArray) ∇iterate
@drule iterate(x::AbstractArray, i::Integer) ∇iterate

function ∇iterate(dy, ::typeof(iterate), t::Tuple)
    return NoTangent(), ungetfield(dy[1], t, 1)
end
function ∇iterate(dy, ::typeof(iterate), t::Tuple, i::Int)
    return NoTangent(), ungetfield(dy[1], t, i), ZeroTangent()
end
@drule iterate(x::Tuple) ∇iterate
@drule iterate(x::Tuple, i::Integer) ∇iterate

# here we explicitely stop propagation in iteration
# over ranges (e.g for i=1:3 ... end)
function ∇iterate(dy, ::typeof(iterate), x::UnitRange)
    return NoTangent(), ZeroTangent()
end
function ∇iterate(dy, ::typeof(iterate), x::UnitRange, i::Int)
    return NoTangent(), ZeroTangent(), ZeroTangent()
end
@drule iterate(x::UnitRange) ∇iterate
@drule iterate(x::UnitRange, i::Integer) ∇iterate

## tuple unpacking

function ∇indexed_iterate(dy, ::typeof(Base.indexed_iterate), t::Tuple, i::Int)
    return NoTangent(), ungetfield(dy[1], t, i), ZeroTangent()
end
@drule Base.indexed_iterate(t::Tuple, i::Int) ∇indexed_iterate

function ∇indexed_iterate(dy, ::typeof(Base.indexed_iterate), t::Tuple, i::Int, state::Int)
    return NoTangent(), ungetfield(dy[1], t, i), ZeroTangent(), ZeroTangent()
end
@drule Base.indexed_iterate(t::Tuple, i::Int, state::Int) ∇indexed_iterate

## tuple construction

∇tuple(dy, ::typeof(tuple), args...) = (NoTangent(), [dy[i] for i=1:length(args)]...)
@drule tuple(args::Vararg) ∇tuple

## some no diff functions

@drule Core.kwfunc(f::Any) (dy, _, f) -> NoTangent()

## cat & co.

function ∇cat_kw(dy, ::typeof(Core.kwfunc(cat)), kw::Any, ::typeof(cat), arrs...)
    return (
        NoTangent(),
        NoTangent(),
        NoTangent(),
        [uncat(dy, i, arrs...; dims=kw.dims) for i=1:length(arrs)]...
    )
end
@drule Core.kwfunc(cat)(kw::Any, _::typeof(cat), arrs::Vararg) ∇cat_kw

function ∇vcat(dy, ::typeof(vcat), arrs...)
    return NoTangent(), [uncat(dy, i, arrs...; dims=1) for i=1:length(arrs)]...
end
@drule vcat(arrs::Vararg) ∇vcat

function ∇hcat(dy, ::typeof(hcat), arrs...)
    return NoTangent(), [uncat(dy, i, arrs...; dims=2) for i=1:length(arrs)]...
end
@drule hcat(arrs::Vararg) ∇hcat


text:
# back propagation
export backprop

function backprop(net::Vector{Layer},inputBase,outputBase;nIter=100)
  if ndims(inputBase)==2
    inputBase=reshape(inputBase,size(inputBase)...,1)
  end
  err=Float64[]
  dataLayerIndex=find(layer->typeof(layer)==JuLearn.dataLayer,net)
  lossLayerIndex=find(layer->typeof(layer)==JuLearn.lossLayer,net)
  nData=size(inputBase)[end-1]
  for iter=1:nIter
    idx=(iter-1)*batchsize+1
    range=map(x->mod1(x,nData),idx:idx+batchsize-1)
    for (i,index) in enumerate(dataLayerIndex)
      net[index].rₒ=inputBase[:,range,i]
      net[index].rᵢ=outputBase[:,range,i]
    end
    forward!(net)
    backward!(net)
    update!(net,iter)
    shareParams!(net)
    push!(err,sum(i->net[i].loss,lossLayerIndex))
    println("iteration=",iter)
  end
   err
end


text:
using Dates

"""
    create_TS

Generates the target time points of the study.
These are the timepoints for which PSCOPF will decide on production levels.
For now, we have 4 taget timepoints starting at a reference timepoint `ts1` and separated by 15 minutes

# Arguments
    - `ts1::Dates.DateTime` : the starting time of the target period of time for the study
"""
function create_target_timepoints(ts1::Dates.DateTime)
    return ts1 .+ [Dates.Minute(0), Dates.Minute(15), Dates.Minute(30), Dates.Minute(45)]
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

abstract type SpacePartitioningAlgorithm end
export SpacePartitioningAlgorithm

"""
	KDTreePartitioning

*BAT-internal, not part of stable public API.*

K-D binary space partitioning algorithm. By default, all parameters are
considered for partitioning. Partition parameters can be specified
manually by using `partition_dims` argument. By default, bounds of the
partitioning tree are extended to those given by prior. This can be changed
by setting `extend_bounds = false`.

Constructor:

	PartitionedSampling(;partition_dims=:auto, extend_bounds::Bool = true)

"""
@with_kw struct KDTreePartitioning <: SpacePartitioningAlgorithm
	partition_dims::Union{Array{Int64,1}, Symbol} = :auto
	extend_bounds::Bool = true
end

"""
	SpacePartTree

*BAT-internal, not part of stable public API.*

The structure stores a partitioning tree generated by any `SpacePartitioningAlgorithm`.
Variables:

* `terminated_leaf` : `true` if the tree node is terminal, `false` otherwise.

* `bounds` : Low and high bound of the tree leaf.

* `left_child` : The left child of the tree, `missing` is the node is terminal.

* `right_child` : The right child of the tree, `missing` is the node is terminal.

* `cut_axis` : Axis along which cut is performed, `missing` if the node is terminal.

* `cut_coordinate` : Coordinate at which a cut is performed, `missing` if the node is terminal.

* `cost` : The sum of the cost functions over leaves.

* `cost_part` : The cost function of the current leaf, `missing` if the node is terminal..
"""
mutable struct SpacePartTree
   terminated_leaf::Bool
   bounds::Array{AbstractFloat}
   left_child::Union{SpacePartTree, Missing}
   right_child::Union{SpacePartTree, Missing}
   cut_axis::Union{Integer, Missing}
   cut_coordinate::Union{AbstractFloat, Missing}
   cost::AbstractFloat
   cost_part::Union{AbstractFloat, Missing}
end


text:
module ASTModels

include("walk1d.jl")
using .Walk1Ds
export
    Walk1DParams,
    Walk1D,
    Walk1DState,
    Walk1DAction

include("gridworld.jl")
using .GridWorlds
export
    GWParams,
    GWSim,
    GWState

end # module


text:
function domain end
function differentiate end
function integrate end
function localdot end

const ∂ = differentiate
const ∫ = integrate

differentiate(a::Union{Number,AbstractArray}) = zero(a)
localdot(a::T, b::T) where {T<:Union{Number,AbstractArray}} = dot(a, b)

struct DomainMismatch <: Exception end

Base.show(io::IO, ::DomainMismatch) =
    Base.print(io, "DomainMismatch(): function space arguments have non-matching domain.")

# FunctionSeries
abstract type FunctionSpace{T} end # functions that support taking linear combinations
abstract type FunctionSeries{T} <: FunctionSpace{T} end
const Const = Union{Number, AbstractArray}

scalartype(::Type{<:FunctionSpace{T}}) where T = scalartype(T)
Base.eltype(::Type{<:FunctionSpace{T}}) where T = T

_rtoldefault(x, y, atol) =
    Base.rtoldefault(scalartype(eltype(x)), scalartype(eltype(y)), atol)

function LinearAlgebra.isapprox(x::FunctionSpace, y::FunctionSpace;
                                atol::Real=0,
                                rtol::Real=_rtoldefault(eltype(x), eltype(y), atol))
    return norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))
end

# AbstractPiecewise
abstract type AbstractPiecewise{T,F<:FunctionSpace{T}} <: FunctionSpace{T} end

domain(p::AbstractPiecewise) = (first(nodes(p)), last(nodes(p)))
function domain(p::AbstractPiecewise, i)
    1 <= i <= length(p) || throw(BoundsError(p, i))
    n = nodes(p)
    @inbounds begin
        return (n[i], n[i+1])
    end
end

Base.eltype(t::AbstractPiecewise{T}) where {T} = eltype(T)
Base.eltype(::Type{<:AbstractPiecewise{T}}) where {T} = eltype(T)


text:
using PyPlot
using WignerFamilies

# compute all j₂=100, j₃=60, m₂=70, m₃=-55, m₁=-m₂-m₃
j₂, j₃, m₂, m₃ = 100, 60, 70, -55
m₁ = -m₂ - m₃
w3j = wigner3j_f(j₂, j₃, m₂, m₃)
js = collect(eachindex(w3j))

plt.clf()
plt.plot(js, w3j.symbols .* 1000)
plt.xlabel(raw"$j$")
plt.ylabel(raw"$f(j) \, \times \, 10^3$")
plt.ylim(-20,20)
plt.axvline(49, ls="--", lw=1)
plt.axvline(98, ls="--", lw=1)
plt.tight_layout()
# plt.savefig("examples/luscombe_and_luban_1998.svg")
gcf()


text:
using SimMdlPrices
using Test

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


text:
export PerfRNN, generate

using StatsBase
using Flux: gate

struct PerfRNN
    model
    encoder::PerformanceOneHotEncoding
    num_velocitybins::Int
end

"""
    generate(perfrnn::PerfRNN, performance::Performance)

Generate `PerformanceEvent`s by sampling from a PerformanceRNN model.
"""
function generate(perfrnn::PerfRNN;
                  primer::Performance=Performance(100, velocity_bins=perfrnn.num_velocitybins),
                  numsteps=3000,
                  as_notesequence=false)

    model, encoder = perfrnn.model, perfrnn.encoder
    Flux.reset!(model)

    performance = deepcopy(primer)

    if isempty(performance)
        push!(performance, encoder.defaultevent)
    end

    # Primer is already numsteps or longer
    if performance.numsteps >= numsteps
        return performance
    end

    indices = map(event -> encode_event(event, encoder), performance)
    inputs = map(index -> Flux.onehot(index, encoder.labels), indices)

    outputs = model.(inputs)
    out = wsample(encoder.labels, softmax(outputs[end]))
    push!(performance, decode_event(out, encoder))

    while performance.numsteps < numsteps
        input = Flux.onehot(out, encoder.labels)
        logits = model(input)
        out = wsample(encoder.labels, softmax(logits))
        push!(performance, decode_event(out, encoder))
    end

    ns = getnotesequence(performance)
    as_notesequence == true && return ns
    midifile(ns)
end

# Replace Flux's LSTMCell with BasicLSTMCell from TensorFlow 1.
# It's implementation is slightly different.
# https://github.com/tensorflow/tensorflow/blob/v2.4.1/tensorflow/python/keras/layers/legacy_rnn/rnn_cell_impl.py#L786-L787
# https://github.com/FluxML/Flux.jl/blob/master/src/layers/recurrent.jl#L141-L142
# The replacement is required because the pre-trained models use BasicLSTMCell.
struct BasicLSTMCell{A,V,S}
  Wi::A
  Wh::A
  b::V
  state0::S
end

function BasicLSTMCell(in::Integer, out::Integer;
                  init = Flux.glorot_uniform,
                  initb = zeros,
                  init_state = zeros)
  cell = BasicLSTMCell(init(out * 4, in), init(out * 4, out), initb(out * 4), (init_state(out,1), init_state(out,1)))
  cell.b[gate(out, 2)] .= 1
  return cell
end

function (m::BasicLSTMCell)((h, c), x) where {A,V,T}
  b, o = m.b, size(h, 1)
  g = m.Wi*x .+ m.Wh*h .+ b
  input = σ.(gate(g, o, 1))
  cell = tanh.(gate(g, o, 2))
  forget = σ.(gate(g, o, 3))
  output = σ.(gate(g, o, 4))
  c = forget .* c .+ input .* cell
  h′ = output .* tanh.(c)
  sz = size(x)
  return (h′, c), reshape(h′, :, sz[2:end]...)
end

Flux.@functor BasicLSTMCell

Base.show(io::IO, l::BasicLSTMCell) =
  print(io, "BasicLSTMCell(", size(l.Wi, 2), ", ", size(l.Wi, 1)÷4, ")")

"""
    BasicLSTM(in::Integer, out::Integer)
Implements the BasicLSTMCell from TensorFlow 1
https://github.com/tensorflow/tensorflow/blob/v2.4.1/tensorflow/python/keras/layers/legacy_rnn/rnn_cell_impl.py#L648-L814
"""
BasicLSTM(a...; ka...) = Flux.Recur(BasicLSTMCell(a...; ka...))
Flux.Recur(m::BasicLSTMCell) = Flux.Recur(m, m.state0)

max_stars_repo_name: VasanthManiVasi/PerformanceRNN

max_issues_repo_name: VasanthManiVasi/PerformanceRNN

max_forks_repo_name: VasanthManiVasi/PerformanceRNN


text:
using PyPlot
plot(1:10)
show()

max_stars_repo_name: szcf-weiya/MCMC

max_issues_repo_name: szcf-weiya/MCMC

max_forks_repo_name: szcf-weiya/MCMC


text:
using PyPlot, JLD2, Polynomials, PyCall, Images
include("../utils/Simulation.jl")
include("../utils/plotting.jl")

n = (128, 128, 128)
JLD2.@load "256Kphi.jld2"

K = K[1:2:end,1:2:end,1:2:end]
phi = phi[1:2:end,1:2:end,1:2:end]

#K = Float32.(imfilter(K, Kernel.gaussian((2,2,2))
d = (50f0, 50f0, 12f0)
qinj = (d[1]*n[1]/2, d[2]*n[2]/2, d[3]*n[3]+(341-256)*6f0-300f0)
qrate = 7
time = 80
nt = 10

sat, p = TwoPhase(K, phi, qinj, qrate, d, time, nt; o=(0f0,0f0,(341-256)*6f0))

JLD2.@save "satpCompass256.jld2" sat p

#figure();imshow(K[:,32,:]')
#figure();imshow(sat[end][:,32,:]')

#plot_3D(reverse(sat[end],dims=3))
#plot_3D(reverse(K,dims=3))
#plot_3D(reverse(phi,dims=3))

max_stars_repo_name: slimgroup/TwoPhaseSpeedTest

max_issues_repo_name: slimgroup/TwoPhaseSpeedTest

max_forks_repo_name: slimgroup/TwoPhaseSpeedTest


text:
using Clang.Generators
using ISL_jll

const ISL_INCLUDE = normpath(joinpath(ISL_jll.artifact_dir, "include"))
const ISL_HEADER_DIR = joinpath(ISL_INCLUDE, "isl")

const BLACKLIST = ["cpp.h", "hmap.h", "hmap_templ.c", "maybe_templ.h", "polynomial_type.h", "val_gmp.h"]
const ISL_HEADERS = [joinpath(ISL_HEADER_DIR, header) for header in readdir(ISL_HEADER_DIR) if header ∉ BLACKLIST]

options = load_options(joinpath(@__DIR__, "wrap.toml"))
args = String["-I$ISL_INCLUDE"]
ctx = create_context(ISL_HEADERS, args, options)

build!(ctx)


text:
function mimestring(mime::Type{<:MIME}, x)
    buf=IOBuffer()
    show(buf, mime(), x)
    String(take!(buf))
end
mimestring(x) = mimestring(MIME"text/plain", x)

struct MultiTerm <: AbstractTerm
    terms::StatsModels.TupleTerm
end
StatsModels.apply_schema(mt::MultiTerm, sch::StatsModels.Schema, Mod::Type) =
    apply_schema.(mt.terms, Ref(sch), Mod)

@testset "terms" begin

    using Statistics

    @testset "concrete_term" begin
        t = term(:aaa)
        ts = term("aaa")
        @test t == ts
        @test string(t) == "aaa"
        @test mimestring(t) == "aaa(unknown)"

        t0 = concrete_term(t, [3, 2, 1])
        @test string(t0) == "aaa"
        @test mimestring(t0) == "aaa(continuous)"
        @test t0.mean == 2.0
        @test t0.var == var([1,2,3])
        @test t0.min == 1.0
        @test t0.max == 3.0

        t1 = concrete_term(t, [:a, :b, :c])
        @test t1.contrasts isa StatsModels.ContrastsMatrix{DummyCoding}
        @test string(t1) == "aaa"
        @test mimestring(t1) == "aaa(DummyCoding:3→2)"

        t3 = concrete_term(t, [:a, :b, :c], DummyCoding())
        @test t3.contrasts isa StatsModels.ContrastsMatrix{DummyCoding}
        @test string(t3) == "aaa"
        @test mimestring(t3) == "aaa(DummyCoding:3→2)"

        t2 = concrete_term(t, [:a, :a, :b], EffectsCoding())
        @test t2.contrasts isa StatsModels.ContrastsMatrix{EffectsCoding}
        @test mimestring(t2) == "aaa(EffectsCoding:2→1)"
        @test string(t2) == "aaa"

        t2full = concrete_term(t, [:a, :a, :b], StatsModels.FullDummyCoding())
        @test t2full.contrasts isa StatsModels.ContrastsMatrix{StatsModels.FullDummyCoding}
        @test mimestring(t2full) == "aaa(StatsModels.FullDummyCoding:2→2)"
        @test string(t2full) == "aaa"
    end

    @testset "term operators" begin
        a = term(:a)
        b = term(:b)
        @test a + b == (a, b)
        @test (a ~ b) == FormulaTerm(a, b)
        @test string(a~b) == "$a ~ $b"
        @test mimestring(a~b) ==
            """FormulaTerm
               Response:
                 a(unknown)
               Predictors:
                 b(unknown)"""
        @test mimestring(a ~ term(1) + b) ==
            """FormulaTerm
               Response:
                 a(unknown)
               Predictors:
                 1
                 b(unknown)"""
        @test a & b == InteractionTerm((a,b))
        @test string(a & b) == "$a & $b"
        @test mimestring(a & b) == "a(unknown) & b(unknown)"
        c = term(:c)
        ab = a+b
        bc = b+c
        abc = a+b+c
        @test ab+c == abc
        @test ab+a == ab
        @test a+bc == abc
        @test b+ab == ab
        @test ab+ab == ab
        @test ab+bc == abc
        @test sum((a,b,c)) == abc
        @test sum((a,)) == a
        @test +a == a
    end

    @testset "uniqueness of FunctionTerms" begin
        f1 = @formula(y ~ lag(x,1) + lag(x,1))
        f2 = @formula(y ~ lag(x,1))
        f3 = @formula(y ~ lag(x,1) + lag(x,2))

        @test f1.rhs == f2.rhs
        @test f1.rhs != f3.rhs

        ## addition of two identical function terms
        @test f2.rhs + f2.rhs == f2.rhs
    end

    @testset "expand nested tuples of terms during apply_schema" begin
        sch = schema((a=rand(10), b=rand(10), c=rand(10)))

        # nested tuples of terms are expanded by apply_schema
        terms = (term(:a), (term(:b), term(:c)))
        terms2 = apply_schema(terms, sch, Nothing)
        @test terms2 isa NTuple{3, ContinuousTerm}
        @test terms2 == apply_schema(term.((:a, :b, :c)), sch, Nothing)

        # a term that generates multiple terms after apply_schema
        mterms = (terms[1], MultiTerm(terms[2]))
        terms3 = apply_schema(mterms, sch, Nothing)

        @test terms2 == terms3
    end

    @testset "Intercept and response traits" begin

        has_responses = [term(:y), term(1), InterceptTerm{true}(), term(:y)+term(:z),
                        term(:y) + term(0), term(:y) + InterceptTerm{false}()]
        no_responses = [term(0), InterceptTerm{false}()]

        has_intercepts = [term(1), InterceptTerm{true}()]
        omits_intercepts = [term(0), term(-1), InterceptTerm{false}()]

        using StatsModels: hasresponse, hasintercept, omitsintercept

        a = term(:a)

        for lhs in has_responses, rhs in has_intercepts
            @test hasresponse(lhs ~ rhs)
            @test hasintercept(lhs ~ rhs)
            @test !omitsintercept(lhs ~ rhs)

            @test hasresponse(lhs ~ rhs + a)
            @test hasintercept(lhs ~ rhs + a)
            @test !omitsintercept(lhs ~ rhs + a)

        end

        for lhs in no_responses, rhs in has_intercepts
            @test !hasresponse(lhs ~ rhs)
            @test hasintercept(lhs ~ rhs)
            @test !omitsintercept(lhs ~ rhs)

            @test !hasresponse(lhs ~ rhs + a)
            @test hasintercept(lhs ~ rhs + a)
            @test !omitsintercept(lhs ~ rhs + a)
        end

        for lhs in has_responses, rhs in omits_intercepts
            @test hasresponse(lhs ~ rhs)
            @test !hasintercept(lhs ~ rhs)
            @test omitsintercept(lhs ~ rhs)

            @test hasresponse(lhs ~ rhs + a)
            @test !hasintercept(lhs ~ rhs + a)
            @test omitsintercept(lhs ~ rhs + a)
        end

        for lhs in no_responses, rhs in omits_intercepts
            @test !hasresponse(lhs ~ rhs)
            @test !hasintercept(lhs ~ rhs)
            @test omitsintercept(lhs ~ rhs)

            @test !hasresponse(lhs ~ rhs + a)
            @test !hasintercept(lhs ~ rhs + a)
            @test omitsintercept(lhs ~ rhs + a)
        end

    end

    @testset "Tuple terms" begin
        using StatsModels: TermOrTerms, TupleTerm, Term
        a, b, c = Term.((:a, :b, :c))

        # TermOrTerms - one or more AbstractTerms (if more, a tuple)
        # empty tuples are never terms
        @test !(() isa TermOrTerms)
        @test (a, ) isa TermOrTerms
        @test (a, b) isa TermOrTerms
        @test (a, b, a&b) isa TermOrTerms
        @test !(((), a) isa TermOrTerms)
        # can't contain further tuples
        @test !((a, (a,), b) isa TermOrTerms)

        # a tuple of AbstractTerms OR Tuples of one or more terms
        # empty tuples are never terms
        @test !(() isa TupleTerm)
        @test (a, ) isa TupleTerm
        @test (a, b) isa TupleTerm
        @test (a, b, a&b) isa TupleTerm
        @test !(((), a) isa TupleTerm)
        @test (((a,), a) isa TupleTerm)

        # no methods for operators on term and empty tuple (=no type piracy)
        @test_throws MethodError a + ()
        @test_throws MethodError () + a
        @test_throws MethodError a & ()
        @test_throws MethodError () & a
        @test_throws MethodError a ~ ()
        @test_throws MethodError () ~ a

        # show methods of empty tuples preserved
        @test "$(())" == "()"
        @test "$((a,b))" == "a + b"
        @test "$((a, ()))" == "(a, ())"
    end

end


text:
abstract type AbstractSquareLattice <:AbstractParallelogram end

export SquareLattice

struct SquareLattice <: AbstractSquareLattice
    Nx::Int
    Ny::Int
end

function Base.getindex(lt::SquareLattice, i::Real, j::Real)
    step = unit(lt)
    (i-0.5)*step, (j-0.5)*step
end

Base.size(sq::SquareLattice) = (sq.Nx, sq.Ny)
vertices(sq::SquareLattice) = 1:sq.Nx*sq.Ny
vec_a(sq::SquareLattice) = [1.0, 0.0]
vec_b(sq::SquareLattice) = [0.0, 1.0]

function isconnected(sq::SquareLattice, i::Int, j::Int)
    u = unit(sq)
    d = distance(sq[i], sq[j])
    d > 0.999u && d < 1.001u
end

function bonds(sq::SquareLattice; nth::Int=1)
    bbb = [1.0, sqrt(2), 2.0, sqrt(5)]
    if nth > length(bbb)
        error("We haven't implemented more than a general n-th order gradient for n>4.
            Please give us a PR or file an issue!")
    end
    d =  bbb[nth] * unit(sq)
    edges = Tuple{Int,Int}[]
    nv = length(sq)
    for i=1:nv
        for j=i+1:nv
            di = distance(sq[i], sq[j])
            if di > d*0.999 && di < d*1.001
                push!(edges, (i, j))
            end
        end
    end
    edges
end

"""
    unit(lattice)

The unit scale of the lattice.
"""
unit(lt::SquareLattice) = 1/(max(lt.Nx, lt.Ny))


text:
function sign!(aws::AbstractAWSConfig, request::Request; time::DateTime=now(Dates.UTC))
    if request.service in ("sdb", "importexport")
        sign_aws2!(aws, request, time)
    else
        sign_aws4!(aws, request, time)
    end
end

function sign_aws2!(aws::AbstractAWSConfig, request::Request, time::DateTime)
    # Create AWS Signature Version 2 Authentication query parameters.
    # http://docs.aws.amazon.com/general/latest/gr/signature-version-2.html

    query = Dict{String,String}()
    for elem in split(request.content, '&'; keepempty=false)
        (n, v) = split(elem, "=")
        query[n] = HTTP.unescapeuri(v)
    end

    request.headers["Content-Type"] = "application/x-www-form-urlencoded; charset=utf-8"

    creds = check_credentials(credentials(aws))
    query["AWSAccessKeyId"] = creds.access_key_id
    query["Expires"] = Dates.format(
        time + Dates.Minute(2), dateformat"yyyy-mm-dd\THH:MM:SS\Z"
    )
    query["SignatureVersion"] = "2"
    query["SignatureMethod"] = "HmacSHA256"

    if !isempty(creds.token)
        query["SecurityToken"] = creds.token
    end

    query = [k => query[k] for k in sort!(collect(keys(query)))]
    uri = HTTP.URI(request.url)
    to_sign = "POST\n$(uri.host)\n$(uri.path)\n$(HTTP.escapeuri(query))"
    push!(
        query,
        "Signature" => strip(base64encode(digest(MD_SHA256, to_sign, creds.secret_key))),
    )

    request.content = HTTP.escapeuri(query)

    return request
end

function sign_aws4!(aws::AbstractAWSConfig, request::Request, time::DateTime)
    # Create AWS Signature Version 4 Authentication Headers.
    # http://docs.aws.amazon.com/general/latest/gr/signature-version-4.html

    date = Dates.format(time, dateformat"yyyymmdd")
    datetime = Dates.format(time, dateformat"yyyymmdd\THHMMSS\Z")

    # Authentication scope...
    authentication_scope = [date, region(aws), request.service, "aws4_request"]

    creds = check_credentials(credentials(aws))
    signing_key = "AWS4$(creds.secret_key)"

    for scope in authentication_scope
        signing_key = digest(MD_SHA256, scope, signing_key)
    end

    # Authentication scope string...
    authentication_scope = join(authentication_scope, "/")

    # SHA256 hash of content...
    content_hash = bytes2hex(digest(MD_SHA256, request.content))

    # HTTP headers...
    delete!(request.headers, "Authorization")

    merge!(
        request.headers,
        Dict(
            "x-amz-content-sha256" => content_hash,
            "x-amz-date" => datetime,
            "Content-MD5" => base64encode(digest(MD_MD5, request.content)),
        ),
    )

    if !isempty(creds.token)
        request.headers["x-amz-security-token"] = creds.token
    end

    # Sort and lowercase() Headers to produce canonical form...
    canonical_headers = join(
        sort!(["$(lowercase(k)):$(strip(v))" for (k, v) in request.headers]), "\n"
    )
    signed_headers = join(sort!([lowercase(k) for k in keys(request.headers)]), ";")

    # Sort Query String...
    uri = HTTP.URI(request.url)
    query = HTTP.URIs.queryparams(uri.query)
    query = [k => query[k] for k in sort!(collect(keys(query)))]

    # Create hash of canonical request...
    canonical_form = string(
        request.request_method,
        "\n",
        request.service == "s3" ? uri.path : HTTP.escapepath(uri.path),
        "\n",
        HTTP.escapeuri(query),
        "\n",
        canonical_headers,
        "\n\n",
        signed_headers,
        "\n",
        content_hash,
    )

    canonical_hash = bytes2hex(digest(MD_SHA256, canonical_form))

    # Create and sign "String to Sign"...
    string_to_sign = "AWS4-HMAC-SHA256\n$datetime\n$authentication_scope\n$canonical_hash"
    signature = bytes2hex(digest(MD_SHA256, string_to_sign, signing_key))

    # Append Authorization header...
    request.headers["Authorization"] = join(
        [
            "AWS4-HMAC-SHA256 Credential=$(creds.access_key_id)/$authentication_scope",
            "SignedHeaders=$signed_headers",
            "Signature=$signature",
        ],
        ", ",
    )

    return request
end


text:
#############################
# probability distributions #
#############################

import Distributions
using SpecialFunctions: lgamma, lbeta, digamma

abstract type Distribution{T} end

"""
    val::T = random(dist::Distribution{T}, args...)

Sample a random choice from the given distribution with the given arguments.
"""
function random end

"""
    lpdf = logpdf(dist::Distribution{T}, value::T, args...)

Evaluate the log probability (density) of the value.
"""
function logpdf end

"""
    has::Bool = has_output_grad(dist::Distribution)

Return true of the gradient if the distribution computes the gradient of the logpdf with respect to the value of the random choice.
"""
function has_output_grad end

"""
    grads::Tuple = logpdf_grad(dist::Distribution{T}, value::T, args...)

Compute the gradient of the logpdf with respect to the value, and each of the arguments.

If `has_output_grad` returns false, then the first element of the returned tuple is `nothing`.
Otherwise, the first element of the tuple is the gradient with respect to the value.
If the return value of `has_argument_grads` has a false value for at position `i`, then the `i+1`th element of the returned tuple has value `nothing`.
Otherwise, this element contains the gradient with respect to the `i`th argument.
"""
function logpdf_grad end

# NOTE: has_argument_grad is documented and exported in gen_fn_interface.jl

get_return_type(::Distribution{T}) where {T} = T

export Distribution
export random
export logpdf
export logpdf_grad
export has_output_grad

include("bernoulli.jl")
include("normal.jl")
include("mvnormal.jl")
include("gamma.jl")
include("inv_gamma.jl")
include("beta.jl")
include("categorical.jl")
include("categorical_with_support.jl")
include("uniform_discrete.jl")
include("uniform_continuous.jl")
include("poisson.jl")
include("piecewise_uniform.jl")
include("beta_uniform.jl")
include("geometric.jl")
include("exponential.jl")

###############
# combinators #
###############

# code shared by vector-shaped combinators
include("vector.jl")

# built-in generative function combinators
include("choice_at/choice_at.jl")
include("call_at/call_at.jl")
include("map/map.jl")
include("unfold/unfold.jl")
include("recurse/recurse.jl")

max_stars_repo_name: rameshputalapattu/Gen

max_issues_repo_name: rameshputalapattu/Gen

max_forks_repo_name: rameshputalapattu/Gen


text:
# Julia program to reverse a number

# Function to reverse a number.
function reverse(num)
    rev = 0
    while(num > 0)
        # Extract the last digit of num and add it to rev
        rev = rev * 10 + num %10
        num = num ÷ 10
    end
    return rev
end

print("Enter the number: ")
num = readline()
num = parse(Int, num)
rev = reverse(num)
println("The reverse of the given number is $rev.")

"""
Time Complexity: O(log(n)), where 'n' is the given number
Space Complexity: O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE 1
Enter the number: 1234
The reverse of the given number is 4321.

SAMPLE 2
Enter the number: 785487
The reverse of the given number is 784587.
"""


text:
using .NumericalFluxes: GradNumericalPenalty, diffusive_boundary_penalty!,
                        diffusive_penalty!,
                        NumericalFluxNonDiffusive, NumericalFluxDiffusive,
                        numerical_flux_nondiffusive!,
                        numerical_boundary_flux_nondiffusive!,
                        numerical_flux_diffusive!,
                        numerical_boundary_flux_diffusive!

using ..Mesh.Geometry

using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

# {{{ FIXME: remove this after we've figure out how to pass through to kernel
const _ξ1x1, _ξ2x1, _ξ3x1 = Grids._ξ1x1, Grids._ξ2x1, Grids._ξ3x1
const _ξ1x2, _ξ2x2, _ξ3x2 = Grids._ξ1x2, Grids._ξ2x2, Grids._ξ3x2
const _ξ1x3, _ξ2x3, _ξ3x3 = Grids._ξ1x3, Grids._ξ2x3, Grids._ξ3x3
const _M, _MI = Grids._M, Grids._MI
const _x1, _x2, _x3 = Grids._x1, Grids._x2, Grids._x3
const _JcV = Grids._JcV

const _n1, _n2, _n3 = Grids._n1, Grids._n2, Grids._n3
const _sM, _vMI = Grids._sM, Grids._vMI
# }}}

"""
    volumerhs!(bl::BalanceLaw, Val(polyorder), rhs, Q, Qvisc, auxstate,
               vgeo, t, D, elems)

Computational kernel: Evaluate the volume integrals on right-hand side of a
`DGBalanceLaw` semi-discretization.

See [`odefun!`](@ref) for usage.
"""
function volumerhs!(bl::BalanceLaw, ::Val{dim}, ::Val{polyorder}, ::direction,
                    rhs, Q, Qvisc, auxstate, vgeo, t,
                    ω, D, elems, increment) where {dim, polyorder, direction}
  N = polyorder
  FT = eltype(Q)
  nstate = num_state(bl,FT)
  nviscstate = num_diffusive(bl,FT)
  nauxstate = num_aux(bl,FT)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  s_F = @shmem FT (3, Nq, Nq, Nqk, nstate)
  s_ω = @shmem FT (Nq, )
  s_half_D = @shmem FT (Nq, Nq)
  l_rhs = @scratch FT (nstate, Nq, Nq, Nqk) 3

  source! !== nothing && (l_S = MArray{Tuple{nstate}, FT}(undef))
  l_Q = @scratch FT (nstate, Nq, Nq, Nqk) 3
  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)
  l_aux = @scratch FT (nauxstate, Nq, Nq, Nqk) 3
  l_F = MArray{Tuple{3, nstate}, FT}(undef)
  l_M = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ1x1 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ1x2 = @scratch FT (Nq, Nq, Nqk) 3
  l_ξ1x3 = @scratch FT (Nq, Nq, Nqk) 3
  if dim == 3 || (dim == 2 && direction == EveryDirection)
    l_ξ2x1 = @scratch FT (Nq, Nq, Nqk) 3
    l_ξ2x2 = @scratch FT (Nq, Nq, Nqk) 3
    l_ξ2x3 = @scratch FT (Nq, Nq, Nqk) 3
  end
  if dim == 3 && direction == EveryDirection
    l_ξ3x1 = @scratch FT (Nq, Nq, Nqk) 3
    l_ξ3x2 = @scratch FT (Nq, Nq, Nqk) 3
    l_ξ3x3 = @scratch FT (Nq, Nq, Nqk) 3
  end

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      s_ω[j] = ω[j]
      @loop for i in (1:Nq; threadIdx().x)
        s_half_D[i, j] = D[i, j] / 2
      end
    end
  end

  @inbounds @views @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          l_M[i, j, k] = vgeo[ijk, _M, e]
          l_ξ1x1[i, j, k] = vgeo[ijk, _ξ1x1, e]
          l_ξ1x2[i, j, k] = vgeo[ijk, _ξ1x2, e]
          l_ξ1x3[i, j, k] = vgeo[ijk, _ξ1x3, e]
          if dim == 3 || (dim == 2 && direction == EveryDirection)
            l_ξ2x1[i, j, k] = vgeo[ijk, _ξ2x1, e]
            l_ξ2x2[i, j, k] = vgeo[ijk, _ξ2x2, e]
            l_ξ2x3[i, j, k] = vgeo[ijk, _ξ2x3, e]
          end
          if dim == 3 && direction == EveryDirection
            l_ξ3x1[i, j, k] = vgeo[ijk, _ξ3x1, e]
            l_ξ3x2[i, j, k] = vgeo[ijk, _ξ3x2, e]
            l_ξ3x3[i, j, k] = vgeo[ijk, _ξ3x3, e]
          end

          @unroll for s = 1:nstate
            l_rhs[s, i, j, k] = increment ? rhs[ijk, s, e] : zero(FT)
          end

          @unroll for s = 1:nstate
            l_Q[s, i, j, k] = Q[ijk, s, e]
          end

          @unroll for s = 1:nauxstate
            l_aux[s, i, j, k] = auxstate[ijk, s, e]
          end

          fill!(l_F, -zero(eltype(l_F)))
          flux_nondiffusive!(bl, Grad{vars_state(bl,FT)}(l_F),
                             Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                             Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)

          @unroll for s = 1:nstate
            s_F[1,i,j,k,s] = l_F[1,s]
            s_F[2,i,j,k,s] = l_F[2,s]
            s_F[3,i,j,k,s] = l_F[3,s]
          end

          # if source! !== nothing
          fill!(l_S, -zero(eltype(l_S)))
          source!(bl, Vars{vars_state(bl,FT)}(l_S),
                  Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                  Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)

          @unroll for s = 1:nstate
            l_rhs[s, i, j, k] += l_S[s]
          end
        end
      end
    end
    @synchronize

    # Weak "outside metrics" derivative
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          @unroll for n = 1:Nq
            @unroll for s = 1:nstate
              Dni = s_half_D[n, i] * s_ω[n] / s_ω[i]
              if dim == 3 || (dim == 2 && direction == EveryDirection)
                Dnj = s_half_D[n, j] * s_ω[n] / s_ω[j]
              end
              if dim == 3 && direction == EveryDirection
                Dnk = s_half_D[n, k] * s_ω[n] / s_ω[k]
              end

              # ξ1-grid lines
              l_rhs[s, i, j, k] += l_ξ1x1[i, j, k] * Dni * s_F[1, n, j, k, s]
              l_rhs[s, i, j, k] += l_ξ1x2[i, j, k] * Dni * s_F[2, n, j, k, s]
              l_rhs[s, i, j, k] += l_ξ1x3[i, j, k] * Dni * s_F[3, n, j, k, s]

              # ξ2-grid lines
              if dim == 3 || (dim == 2 && direction == EveryDirection)
                l_rhs[s, i, j, k] += l_ξ2x1[i, j, k] * Dnj * s_F[1, i, n, k, s]
                l_rhs[s, i, j, k] += l_ξ2x2[i, j, k] * Dnj * s_F[2, i, n, k, s]
                l_rhs[s, i, j, k] += l_ξ2x3[i, j, k] * Dnj * s_F[3, i, n, k, s]
              end

              # ξ3-grid lines
              if dim == 3 && direction == EveryDirection
                l_rhs[s, i, j, k] += l_ξ3x1[i, j, k] * Dnk * s_F[1, i, j, n, s]
                l_rhs[s, i, j, k] += l_ξ3x2[i, j, k] * Dnk * s_F[2, i, j, n, s]
                l_rhs[s, i, j, k] += l_ξ3x3[i, j, k] * Dnk * s_F[3, i, j, n, s]
              end
            end
          end
        end
      end
    end
    @synchronize

    # Add in the diffusive flux (multiply by 2 since derivative is halfed)
    # This allows symmetric treament of the 2nd order derivative terms
    # as well as build "inside metrics" flux
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))

          @unroll for s = 1:nviscstate
            l_Qvisc[s] = Qvisc[ijk, s, e]
          end

          fill!(l_F, -zero(eltype(l_F)))
          flux_diffusive!(bl, Grad{vars_state(bl,FT)}(l_F),
                          Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                          Vars{vars_diffusive(bl,FT)}(l_Qvisc),
                          Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)

          @unroll for s = 1:nstate
            F1, F2, F3 = s_F[1,i,j,k,s], s_F[2,i,j,k,s], s_F[3,i,j,k,s]

            F1 += 2l_F[1,s]
            F2 += 2l_F[2,s]
            F3 += 2l_F[3,s]

            s_F[1,i,j,k,s] = l_M[i, j, k] * (l_ξ1x1[i, j, k] * F1 +
                                              l_ξ1x2[i, j, k] * F2 +
                                              l_ξ1x3[i, j, k] * F3)
            if dim == 3 || (dim == 2 && direction == EveryDirection)
              s_F[2,i,j,k,s] = l_M[i, j, k] * (l_ξ2x1[i, j, k] * F1 +
                                               l_ξ2x2[i, j, k] * F2 +
                                               l_ξ2x3[i, j, k] * F3)
            end
            if dim == 3 && direction == EveryDirection
              s_F[3,i,j,k,s] = l_M[i, j, k] * (l_ξ3x1[i, j, k] * F1 +
                                               l_ξ3x2[i, j, k] * F2 +
                                               l_ξ3x3[i, j, k] * F3)
            end
          end
        end
      end
    end
    @synchronize

    # Weak "inside metrics" derivative
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          MI = vgeo[ijk, _MI, e]
          @unroll for s = 1:nstate
            @unroll for n = 1:Nq
              # ξ1-grid lines
              l_rhs[s, i, j, k] += MI * s_half_D[n, i] * s_F[1, n, j, k, s]

              # ξ2-grid lines
              if dim == 3 || (dim == 2 && direction == EveryDirection)
                l_rhs[s, i, j, k] += MI * s_half_D[n, j] * s_F[2, i, n, k, s]
              end

              # ξ3-grid lines
              if dim == 3 && direction == EveryDirection
                l_rhs[s, i, j, k] += MI * s_half_D[n, k] * s_F[3, i, j, n, s]
              end
            end
          end
        end
      end
    end
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          @unroll for s = 1:nstate
            rhs[ijk, s, e] = l_rhs[s, i, j, k]
          end
        end
      end
    end
    @synchronize
  end
  nothing
end

function volumerhs!(bl::BalanceLaw, ::Val{dim}, ::Val{polyorder},
                    ::VerticalDirection, rhs, Q, Qvisc, auxstate, vgeo, t,
                    ω, D, elems, increment) where {dim, polyorder}
  N = polyorder
  FT = eltype(Q)
  nstate = num_state(bl,FT)
  nviscstate = num_diffusive(bl,FT)
  nauxstate = num_aux(bl,FT)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  s_F = @shmem FT (3, Nq, Nq, Nqk, nstate)
  s_ω = @shmem FT (Nq, )
  s_half_D = @shmem FT (Nq, Nq)
  l_rhs = @scratch FT (nstate, Nq, Nq, Nqk) 3

  source! !== nothing && (l_S = MArray{Tuple{nstate}, FT}(undef))
  l_Q = @scratch FT (nstate, Nq, Nq, Nqk) 3
  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)
  l_aux = @scratch FT (nauxstate, Nq, Nq, Nqk) 3
  l_F = MArray{Tuple{3, nstate}, FT}(undef)
  l_M = @scratch FT (Nq, Nq, Nqk) 3

  l_ζx1 = @scratch FT (Nq, Nq, Nqk) 3
  l_ζx2 = @scratch FT (Nq, Nq, Nqk) 3
  l_ζx3 = @scratch FT (Nq, Nq, Nqk) 3

  _ζx1 = dim == 2 ? _ξ2x1 : _ξ3x1
  _ζx2 = dim == 2 ? _ξ2x2 : _ξ3x2
  _ζx3 = dim == 2 ? _ξ2x3 : _ξ3x3

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      s_ω[j] = ω[j]
      @loop for i in (1:Nq; threadIdx().x)
        s_half_D[i, j] = D[i, j] / 2
      end
    end
  end

  @inbounds @views @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          l_M[i, j, k] = vgeo[ijk, _M, e]
          l_ζx1[i, j, k] = vgeo[ijk, _ζx1, e]
          l_ζx2[i, j, k] = vgeo[ijk, _ζx2, e]
          l_ζx3[i, j, k] = vgeo[ijk, _ζx3, e]

          @unroll for s = 1:nstate
            l_rhs[s, i, j, k] = increment ? rhs[ijk, s, e] : zero(FT)
          end

          @unroll for s = 1:nstate
            l_Q[s, i, j, k] = Q[ijk, s, e]
          end

          @unroll for s = 1:nauxstate
            l_aux[s, i, j, k] = auxstate[ijk, s, e]
          end

          fill!(l_F, -zero(eltype(l_F)))
          flux_nondiffusive!(bl, Grad{vars_state(bl,FT)}(l_F),
                             Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                             Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)

          @unroll for s = 1:nstate
            s_F[1,i,j,k,s] = l_F[1,s]
            s_F[2,i,j,k,s] = l_F[2,s]
            s_F[3,i,j,k,s] = l_F[3,s]
          end

          # if source! !== nothing
          fill!(l_S, -zero(eltype(l_S)))
          source!(bl, Vars{vars_state(bl,FT)}(l_S),
                  Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                  Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)

          @unroll for s = 1:nstate
            l_rhs[s, i, j, k] += l_S[s]
          end
        end
      end
    end
    @synchronize

    # Weak "outside metrics" derivative
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          @unroll for n = 1:Nq
            @unroll for s = 1:nstate
              if dim == 2
                Dnj = s_half_D[n, j] * s_ω[n] / s_ω[j]
                l_rhs[s, i, j, k] += l_ζx1[i, j, k] * Dnj * s_F[1, i, n, k, s]
                l_rhs[s, i, j, k] += l_ζx2[i, j, k] * Dnj * s_F[2, i, n, k, s]
                l_rhs[s, i, j, k] += l_ζx3[i, j, k] * Dnj * s_F[3, i, n, k, s]
              else
                Dnk = s_half_D[n, k] * s_ω[n] / s_ω[k]
                l_rhs[s, i, j, k] += l_ζx1[i, j, k] * Dnk * s_F[1, i, j, n, s]
                l_rhs[s, i, j, k] += l_ζx2[i, j, k] * Dnk * s_F[2, i, j, n, s]
                l_rhs[s, i, j, k] += l_ζx3[i, j, k] * Dnk * s_F[3, i, j, n, s]
              end
            end
          end
        end
      end
    end
    @synchronize

    # Build "inside metrics" flux
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))

          @unroll for s = 1:nviscstate
            l_Qvisc[s] = Qvisc[ijk, s, e]
          end

          fill!(l_F, -zero(eltype(l_F)))
          flux_diffusive!(bl, Grad{vars_state(bl,FT)}(l_F),
                          Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                          Vars{vars_diffusive(bl,FT)}(l_Qvisc),
                          Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)

          @unroll for s = 1:nstate
            F1, F2, F3 = s_F[1,i,j,k,s], s_F[2,i,j,k,s], s_F[3,i,j,k,s]

            F1 += 2l_F[1,s]
            F2 += 2l_F[2,s]
            F3 += 2l_F[3,s]

            s_F[3,i,j,k,s] = l_M[i, j, k] * (l_ζx1[i, j, k] * F1 +
                                             l_ζx2[i, j, k] * F2 +
                                             l_ζx3[i, j, k] * F3)
          end
        end
      end
    end
    @synchronize

    # Weak "inside metrics" derivative
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          MI = vgeo[ijk, _MI, e]
          @unroll for s = 1:nstate
            @unroll for n = 1:Nq
              if dim == 2
                Dnj = s_half_D[n, j]
                l_rhs[s, i, j, k] += MI * Dnj * s_F[3, i, n, k, s]
              else
                Dnk = s_half_D[n, k]
                l_rhs[s, i, j, k] += MI * Dnk * s_F[3, i, j, n, s]
              end
            end
          end
        end
      end
    end
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          @unroll for s = 1:nstate
            rhs[ijk, s, e] = l_rhs[s, i, j, k]
          end
        end
      end
    end
    @synchronize
  end
  nothing
end

"""
    facerhs!(bl::BalanceLaw, Val(polyorder),
            numfluxnondiff::NumericalFluxNonDiffusive,
            numfluxdiff::NumericalFluxDiffusive,
            rhs, Q, Qvisc, auxstate,
            vgeo, sgeo, t, vmapM, vmapP, elemtobndy,
            elems)

Computational kernel: Evaluate the surface integrals on right-hand side of a
`BalanceLaw` semi-discretization.

See [`odefun!`](@ref) for usage.
"""
function facerhs!(bl::BalanceLaw, ::Val{dim}, ::Val{polyorder}, ::direction,
                  numfluxnondiff::NumericalFluxNonDiffusive,
                  numfluxdiff::NumericalFluxDiffusive,
                  rhs, Q, Qvisc, auxstate, vgeo, sgeo, t, vmapM, vmapP,
                  elemtobndy, elems) where {dim, polyorder, direction}
  N = polyorder
  FT = eltype(Q)
  nstate = num_state(bl,FT)
  nviscstate = num_diffusive(bl,FT)
  nauxstate = num_aux(bl,FT)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

  faces = 1:nface
  if direction == VerticalDirection
    faces = nface-1:nface
  elseif direction == HorizontalDirection
    faces = 1:nface-2
  end

  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq

  l_QM = MArray{Tuple{nstate}, FT}(undef)
  l_QviscM = MArray{Tuple{nviscstate}, FT}(undef)
  l_auxM = MArray{Tuple{nauxstate}, FT}(undef)

  # Need two copies since numerical_flux_nondiffusive! can modify QP
  l_QPnondiff = MArray{Tuple{nstate}, FT}(undef)
  l_QPdiff = MArray{Tuple{nstate}, FT}(undef)

  # Need two copies since numerical_flux_nondiffusive! can modify auxP
  l_auxPnondiff = MArray{Tuple{nauxstate}, FT}(undef)
  l_auxPdiff = MArray{Tuple{nauxstate}, FT}(undef)

  l_QviscP = MArray{Tuple{nviscstate}, FT}(undef)

  l_Q_bot1 = MArray{Tuple{nstate}, FT}(undef)
  l_Qvisc_bot1 = MArray{Tuple{nviscstate}, FT}(undef)
  l_aux_bot1 = MArray{Tuple{nauxstate}, FT}(undef)

  l_F = MArray{Tuple{nstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    for f in faces
      @loop for n in (1:Nfp; threadIdx().x)
        nM = SVector(sgeo[_n1, n, f, e], sgeo[_n2, n, f, e], sgeo[_n3, n, f, e])
        sM, vMI = sgeo[_sM, n, f, e], sgeo[_vMI, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        # Load minus side data
        @unroll for s = 1:nstate
          l_QM[s] = Q[vidM, s, eM]
        end

        @unroll for s = 1:nviscstate
          l_QviscM[s] = Qvisc[vidM, s, eM]
        end

        @unroll for s = 1:nauxstate
          l_auxM[s] = auxstate[vidM, s, eM]
        end

        # Load plus side data
        @unroll for s = 1:nstate
          l_QPdiff[s] = l_QPnondiff[s] = l_QPdiff[s] = Q[vidP, s, eP]
        end

        @unroll for s = 1:nviscstate
          l_QviscP[s] = Qvisc[vidP, s, eP]
        end

        @unroll for s = 1:nauxstate
          l_auxPdiff[s] = l_auxPnondiff[s] = auxstate[vidP, s, eP]
        end

        bctype = elemtobndy[f, e]
        fill!(l_F, -zero(eltype(l_F)))
        if bctype == 0
          numerical_flux_nondiffusive!(numfluxnondiff, bl, l_F, nM, l_QM,
                                       l_auxM, l_QPnondiff, l_auxPnondiff, t)
          numerical_flux_diffusive!(numfluxdiff, bl, l_F, nM, l_QM, l_QviscM,
                                    l_auxM, l_QPdiff, l_QviscP, l_auxPdiff, t)
        else
          if (dim == 2 && f == 3) || (dim == 3 && f == 5)
            # Loop up the first element along all horizontal elements
            @unroll for s = 1:nstate
              l_Q_bot1[s] = Q[n + Nqk^2, s, e]
            end
            @unroll for s = 1:nviscstate
              l_Qvisc_bot1[s] = Qvisc[n + Nqk^2, s, e]
            end
            @unroll for s = 1:nauxstate
              l_aux_bot1[s] = auxstate[n + Nqk^2,s, e]
            end
          end
          numerical_boundary_flux_nondiffusive!(numfluxnondiff, bl, l_F, nM,
                                                l_QM, l_auxM, l_QPnondiff,
                                                l_auxPnondiff, bctype, t,
                                                l_Q_bot1, l_aux_bot1)
          numerical_boundary_flux_diffusive!(numfluxdiff, bl, l_F, nM, l_QM,
                                             l_QviscM, l_auxM, l_QPdiff,
                                             l_QviscP, l_auxPdiff, bctype, t,
                                             l_Q_bot1, l_Qvisc_bot1, l_aux_bot1)
        end

        #Update RHS
        @unroll for s = 1:nstate
          # FIXME: Should we pretch these?
          rhs[vidM, s, eM] -= vMI * sM * l_F[s]
        end
      end
      # Need to wait after even faces to avoid race conditions
      f % 2 == 0 && @synchronize
    end
  end
  nothing
end

function volumeviscterms!(bl::BalanceLaw, ::Val{dim}, ::Val{polyorder},
                          ::direction, Q, Qvisc, auxstate, vgeo, t, D,
                          elems) where {dim, polyorder, direction}
  N = polyorder

  FT = eltype(Q)
  nstate = num_state(bl,FT)
  ngradstate = num_gradient(bl,FT)
  nviscstate = num_diffusive(bl,FT)
  nauxstate = num_aux(bl,FT)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  ngradtransformstate = nstate

  s_G = @shmem FT (Nq, Nq, Nqk, ngradstate)
  s_D = @shmem FT (Nq, Nq)

  l_Q = @scratch FT (ngradtransformstate, Nq, Nq, Nqk) 3
  l_aux = @scratch FT (nauxstate, Nq, Nq, Nqk) 3
  l_G = MArray{Tuple{ngradstate}, FT}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)
  l_gradG = MArray{Tuple{3, ngradstate}, FT}(undef)

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        s_D[i, j] = D[i, j]
      end
    end
  end

  @inbounds @views @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          @unroll for s = 1:ngradtransformstate
            l_Q[s, i, j, k] = Q[ijk, s, e]
          end

          @unroll for s = 1:nauxstate
            l_aux[s, i, j, k] = auxstate[ijk, s, e]
          end

          fill!(l_G, -zero(eltype(l_G)))
          gradvariables!(bl, Vars{vars_gradient(bl,FT)}(l_G), Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                     Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)
          @unroll for s = 1:ngradstate
            s_G[i, j, k, s] = l_G[s]
          end
        end
      end
    end
    @synchronize

    # Compute gradient of each state
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          ξ1x1, ξ1x2, ξ1x3 = vgeo[ijk, _ξ1x1, e], vgeo[ijk, _ξ1x2, e], vgeo[ijk, _ξ1x3, e]
          if dim == 3 || (dim == 2 && direction == EveryDirection)
            ξ2x1, ξ2x2, ξ2x3 = vgeo[ijk, _ξ2x1, e], vgeo[ijk, _ξ2x2, e], vgeo[ijk, _ξ2x3, e]
          end
          if dim == 3 && direction == EveryDirection
            ξ3x1, ξ3x2, ξ3x3 = vgeo[ijk, _ξ3x1, e], vgeo[ijk, _ξ3x2, e], vgeo[ijk, _ξ3x3, e]
          end

          @unroll for s = 1:ngradstate
            Gξ1 = Gξ2 = Gξ3 = zero(FT)
            @unroll for n = 1:Nq
              Gξ1 += s_D[i, n] * s_G[n, j, k, s]
              if dim == 3 || (dim == 2 && direction == EveryDirection)
                Gξ2 += s_D[j, n] * s_G[i, n, k, s]
              end
              if dim == 3 && direction == EveryDirection
                Gξ3 += s_D[k, n] * s_G[i, j, n, s]
              end
            end
            l_gradG[1, s] = ξ1x1 * Gξ1
            l_gradG[2, s] = ξ1x2 * Gξ1
            l_gradG[3, s] = ξ1x3 * Gξ1

            if dim == 3 || (dim == 2 && direction == EveryDirection)
              l_gradG[1, s] += ξ2x1 * Gξ2
              l_gradG[2, s] += ξ2x2 * Gξ2
              l_gradG[3, s] += ξ2x3 * Gξ2
            end

            if dim == 3 && direction == EveryDirection
              l_gradG[1, s] += ξ3x1 * Gξ3
              l_gradG[2, s] += ξ3x2 * Gξ3
              l_gradG[3, s] += ξ3x3 * Gξ3
            end
          end

          fill!(l_Qvisc, -zero(eltype(l_Qvisc)))
          diffusive!(bl, Vars{vars_diffusive(bl,FT)}(l_Qvisc), Grad{vars_gradient(bl,FT)}(l_gradG),
                     Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]), Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)

          @unroll for s = 1:nviscstate
            Qvisc[ijk, s, e] = l_Qvisc[s]
          end
        end
      end
    end
    @synchronize
  end
end

function volumeviscterms!(bl::BalanceLaw, ::Val{dim}, ::Val{polyorder},
                          ::VerticalDirection, Q, Qvisc, auxstate, vgeo, t, D,
                          elems) where {dim, polyorder}
  N = polyorder

  FT = eltype(Q)
  nstate = num_state(bl,FT)
  ngradstate = num_gradient(bl,FT)
  nviscstate = num_diffusive(bl,FT)
  nauxstate = num_aux(bl,FT)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  ngradtransformstate = nstate

  s_G = @shmem FT (Nq, Nq, Nqk, ngradstate)
  s_D = @shmem FT (Nq, Nq)

  l_Q = @scratch FT (ngradtransformstate, Nq, Nq, Nqk) 3
  l_aux = @scratch FT (nauxstate, Nq, Nq, Nqk) 3
  l_G = MArray{Tuple{ngradstate}, FT}(undef)
  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)
  l_gradG = MArray{Tuple{3, ngradstate}, FT}(undef)

  _ζx1 = dim == 2 ? _ξ2x1 : _ξ3x1
  _ζx2 = dim == 2 ? _ξ2x2 : _ξ3x2
  _ζx3 = dim == 2 ? _ξ2x3 : _ξ3x3

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for j in (1:Nq; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        s_D[i, j] = D[i, j]
      end
    end
  end

  @inbounds @views @loop for e in (elems; blockIdx().x)
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          @unroll for s = 1:ngradtransformstate
            l_Q[s, i, j, k] = Q[ijk, s, e]
          end

          @unroll for s = 1:nauxstate
            l_aux[s, i, j, k] = auxstate[ijk, s, e]
          end

          fill!(l_G, -zero(eltype(l_G)))
          gradvariables!(bl, Vars{vars_gradient(bl,FT)}(l_G),
                         Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                         Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)
          @unroll for s = 1:ngradstate
            s_G[i, j, k, s] = l_G[s]
          end
        end
      end
    end
    @synchronize

    # Compute gradient of each state
    @loop for k in (1:Nqk; threadIdx().z)
      @loop for j in (1:Nq; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          ijk = i + Nq * ((j-1) + Nq * (k-1))
          ζx1 = vgeo[ijk, _ζx1, e]
          ζx2 = vgeo[ijk, _ζx2, e]
          ζx3 = vgeo[ijk, _ζx3, e]

          @unroll for s = 1:ngradstate
            Gζ = zero(FT)
            @unroll for n = 1:Nq
              if dim == 2
                Gζ += s_D[j, n] * s_G[i, n, k, s]
              elseif dim == 3
                Gζ += s_D[k, n] * s_G[i, j, n, s]
              end
            end
            l_gradG[1, s] = ζx1 * Gζ
            l_gradG[2, s] = ζx2 * Gζ
            l_gradG[3, s] = ζx3 * Gζ
          end

          fill!(l_Qvisc, -zero(eltype(l_Qvisc)))
          diffusive!(bl, Vars{vars_diffusive(bl,FT)}(l_Qvisc),
                     Grad{vars_gradient(bl,FT)}(l_gradG),
                     Vars{vars_state(bl,FT)}(l_Q[:, i, j, k]),
                     Vars{vars_aux(bl,FT)}(l_aux[:, i, j, k]), t)

          @unroll for s = 1:nviscstate
            Qvisc[ijk, s, e] = l_Qvisc[s]
          end
        end
      end
    end
    @synchronize
  end
end

function faceviscterms!(bl::BalanceLaw, ::Val{dim}, ::Val{polyorder},
                        ::direction, gradnumpenalty::GradNumericalPenalty, Q,
                        Qvisc, auxstate, vgeo, sgeo, t, vmapM, vmapP,
                        elemtobndy, elems) where {dim, polyorder, direction}
  N = polyorder
  FT = eltype(Q)
  nstate = num_state(bl,FT)
  ngradstate = num_gradient(bl,FT)
  nviscstate = num_diffusive(bl,FT)
  nauxstate = num_aux(bl,FT)

  if dim == 1
    Np = (N+1)
    Nfp = 1
    nface = 2
  elseif dim == 2
    Np = (N+1) * (N+1)
    Nfp = (N+1)
    nface = 4
  elseif dim == 3
    Np = (N+1) * (N+1) * (N+1)
    Nfp = (N+1) * (N+1)
    nface = 6
  end

  faces = 1:nface
  if direction == VerticalDirection
    faces = nface-1:nface
  elseif direction == HorizontalDirection
    faces = 1:nface-2
  end

  Nqk = dim == 2 ? 1 : N+1

  ngradtransformstate = nstate

  l_QM = MArray{Tuple{ngradtransformstate}, FT}(undef)
  l_auxM = MArray{Tuple{nauxstate}, FT}(undef)
  l_GM = MArray{Tuple{ngradstate}, FT}(undef)

  l_QP = MArray{Tuple{ngradtransformstate}, FT}(undef)
  l_auxP = MArray{Tuple{nauxstate}, FT}(undef)
  l_GP = MArray{Tuple{ngradstate}, FT}(undef)

  l_Qvisc = MArray{Tuple{nviscstate}, FT}(undef)

  l_Q_bot1 = MArray{Tuple{nstate}, FT}(undef)
  l_aux_bot1 = MArray{Tuple{nauxstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    for f in faces
      @loop for n in (1:Nfp; threadIdx().x)
        nM = SVector(sgeo[_n1, n, f, e], sgeo[_n2, n, f, e], sgeo[_n3, n, f, e])
        sM, vMI = sgeo[_sM, n, f, e], sgeo[_vMI, n, f, e]
        idM, idP = vmapM[n, f, e], vmapP[n, f, e]

        eM, eP = e, ((idP - 1) ÷ Np) + 1
        vidM, vidP = ((idM - 1) % Np) + 1,  ((idP - 1) % Np) + 1

        # Load minus side data
        @unroll for s = 1:ngradtransformstate
          l_QM[s] = Q[vidM, s, eM]
        end

        @unroll for s = 1:nauxstate
          l_auxM[s] = auxstate[vidM, s, eM]
        end

        fill!(l_GM, -zero(eltype(l_GM)))
        gradvariables!(bl, Vars{vars_gradient(bl,FT)}(l_GM),
                       Vars{vars_state(bl,FT)}(l_QM),
                       Vars{vars_aux(bl,FT)}(l_auxM), t)

        # Load plus side data
        @unroll for s = 1:ngradtransformstate
          l_QP[s] = Q[vidP, s, eP]
        end

        @unroll for s = 1:nauxstate
          l_auxP[s] = auxstate[vidP, s, eP]
        end

        fill!(l_GP, -zero(eltype(l_GP)))
        gradvariables!(bl, Vars{vars_gradient(bl,FT)}(l_GP),
                       Vars{vars_state(bl,FT)}(l_QP),
                       Vars{vars_aux(bl,FT)}(l_auxP), t)

        bctype = elemtobndy[f, e]
        fill!(l_Qvisc, -zero(eltype(l_Qvisc)))
        if bctype == 0
          diffusive_penalty!(gradnumpenalty, bl, l_Qvisc, nM, l_GM, l_QM,
                             l_auxM, l_GP, l_QP, l_auxP, t)
        else
          if (dim == 2 && f == 3) || (dim == 3 && f == 5)
            # Loop up the first element along all horizontal elements
            @unroll for s = 1:nstate
              l_Q_bot1[s] = Q[n + Nqk^2, s, e]
            end
            @unroll for s = 1:nauxstate
              l_aux_bot1[s] = auxstate[n + Nqk^2,s, e]
            end
          end
          diffusive_boundary_penalty!(gradnumpenalty, bl, l_Qvisc, nM, l_GM,
                                      l_QM, l_auxM, l_GP, l_QP, l_auxP, bctype,
                                      t, l_Q_bot1, l_aux_bot1)
        end

        @unroll for s = 1:nviscstate
          Qvisc[vidM, s, eM] += vMI * sM * l_Qvisc[s]
        end
      end
      # Need to wait after even faces to avoid race conditions
      f % 2 == 0 && @synchronize
    end
  end
  nothing
end

function initstate!(bl::BalanceLaw, ::Val{dim}, ::Val{polyorder}, state, auxstate, vgeo, elems, args...) where {dim, polyorder}
  N = polyorder
  FT = eltype(auxstate)
  nauxstate = num_aux(bl,FT)
  nstate = num_state(bl,FT)

  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq
  Np = Nq * Nq * Nqk

  l_state = MArray{Tuple{nstate}, FT}(undef)
  l_aux = MArray{Tuple{nauxstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for n in (1:Np; threadIdx().x)
      coords = SVector(vgeo[n, _x1, e], vgeo[n, _x2, e], vgeo[n, _x3, e])
      @unroll for s = 1:nauxstate
        l_aux[s] = auxstate[n, s, e]
      end
      @unroll for s = 1:nstate
        l_state[s] = state[n, s, e]
      end
      init_state!(bl, Vars{vars_state(bl,FT)}(l_state), Vars{vars_aux(bl,FT)}(l_aux), coords, args...)
      @unroll for s = 1:nstate
        state[n, s, e] = l_state[s]
      end
    end
  end
end

"""
    initauxstate!(bl::BalanceLaw, Val(polyorder), auxstate, vgeo, elems)

Computational kernel: Initialize the auxiliary state

See [`DGBalanceLaw`](@ref) for usage.
"""
function initauxstate!(bl::BalanceLaw, ::Val{dim}, ::Val{polyorder}, auxstate, vgeo, elems) where {dim, polyorder}
  N = polyorder
  FT = eltype(auxstate)
  nauxstate = num_aux(bl,FT)

  Nq = N + 1
  Nqk = dim == 2 ? 1 : Nq
  Np = Nq * Nq * Nqk

  l_aux = MArray{Tuple{nauxstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for n in (1:Np; threadIdx().x)
      @unroll for s = 1:nauxstate
        l_aux[s] = auxstate[n, s, e]
      end

      init_aux!(bl, Vars{vars_aux(bl,FT)}(l_aux), LocalGeometry(Val(polyorder),vgeo,n,e))

      @unroll for s = 1:nauxstate
        auxstate[n, s, e] = l_aux[s]
      end
    end
  end
end

"""
    knl_nodal_update_aux!(bl::BalanceLaw, ::Val{dim}, ::Val{N}, f!, Q, auxstate,
                          t, elems) where {dim, N}

Update the auxiliary state array
"""
function knl_nodal_update_aux!(bl::BalanceLaw, ::Val{dim}, ::Val{N}, f!, Q,
                               auxstate, t, elems) where {dim, N}
  FT = eltype(Q)
  nstate = num_state(bl,FT)
  nviscstate = num_diffusive(bl,FT)
  nauxstate = num_aux(bl,FT)

  Nq = N + 1

  Nqk = dim == 2 ? 1 : Nq

  Np = Nq * Nq * Nqk

  l_Q = MArray{Tuple{nstate}, FT}(undef)
  l_aux = MArray{Tuple{nauxstate}, FT}(undef)

  @inbounds @loop for e in (elems; blockIdx().x)
    @loop for n in (1:Np; threadIdx().x)
      @unroll for s = 1:nstate
        l_Q[s] = Q[n, s, e]
      end

      @unroll for s = 1:nauxstate
        l_aux[s] = auxstate[n, s, e]
      end

      f!(bl, Vars{vars_state(bl,FT)}(l_Q),
         Vars{vars_aux(bl,FT)}(l_aux), t)

      @unroll for s = 1:nauxstate
        auxstate[n, s, e] = l_aux[s]
      end
    end
  end
end

"""
    knl_indefinite_stack_integral!(::Val{dim}, ::Val{N}, ::Val{nstate},
                                            ::Val{nauxstate}, ::Val{nvertelem},
                                            int_knl!, Q, auxstate, vgeo, Imat,
                                            elems, ::Val{outstate}
                                           ) where {dim, N, nstate, nauxstate,
                                                    outstate, nvertelem}

Computational kernel: compute indefinite integral along the vertical stack

See [`DGBalanceLaw`](@ref) for usage.
"""
function knl_indefinite_stack_integral!(bl::BalanceLaw, ::Val{dim}, ::Val{N}, ::Val{nvertelem},
                                        Q, auxstate, vgeo, Imat,
                                        elems, ::Val{nout}
                                       ) where {dim, N, nvertelem,
                                                nout}
  FT = eltype(Q)
  nstate = num_state(bl,FT)
  nauxstate = num_aux(bl,FT)

  Nq = N + 1
  Nqj = dim == 2 ? 1 : Nq

  l_Q = MArray{Tuple{nstate}, FT}(undef)
  l_aux = MArray{Tuple{nauxstate}, FT}(undef)
  l_knl = MArray{Tuple{nout, Nq}, FT}(undef)
  # note that k is the second not 4th index (since this is scratch memory and k
  # needs to be persistent across threads)
  l_int = @scratch FT (nout, Nq, Nq, Nqj) 2

  s_I = @shmem FT (Nq, Nq)

  @inbounds @loop for k in (1; threadIdx().z)
    @loop for i in (1:Nq; threadIdx().x)
      @unroll for n = 1:Nq
        s_I[i, n] = Imat[i, n]
      end
    end
  end
  @synchronize

  @inbounds @loop for eh in (elems; blockIdx().x)
    # Initialize the constant state at zero
    @loop for j in (1:Nqj; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        @unroll for k in 1:Nq
          @unroll for s = 1:nout
            l_int[s, k, i, j] = 0
          end
        end
      end
    end
    # Loop up the stack of elements
    for ev = 1:nvertelem
      e = ev + (eh - 1) * nvertelem

      # Evaluate the integral kernel at each DOF in the slabk
      @loop for j in (1:Nqj; threadIdx().y)
        @loop for i in (1:Nq; threadIdx().x)
          # loop up the pencil
          @unroll for k in 1:Nq
            ijk = i + Nq * ((j-1) + Nqj * (k-1))
            Jc = vgeo[ijk, _JcV, e]
            @unroll for s = 1:nstate
              l_Q[s] = Q[ijk, s, e]
            end

            @unroll for s = 1:nauxstate
              l_aux[s] = auxstate[ijk, s, e]
            end

            integrate_aux!(bl, Vars{vars_integrals(bl, FT)}(view(l_knl, :, k)),
              Vars{vars_state(bl, FT)}(l_Q), Vars{vars_aux(bl,FT)}(l_aux))

            # multiply in the curve jacobian
            @unroll for s = 1:nout
              l_knl[s, k] *= Jc
            end
          end

          # Evaluate the integral up the element
          @unroll for s = 1:nout
            @unroll for k in 1:Nq
              @unroll for n in 1:Nq
                l_int[s, k, i, j] += s_I[k, n] * l_knl[s, n]
              end
            end
          end

          # Store out to memory and reset the background value for next element
          @unroll for k in 1:Nq
            ijk = i + Nq * ((j-1) + Nqj * (k-1))
            @unroll for ind_out = 1:nout
              auxstate[ijk, ind_out, e] = l_int[ind_out, k, i, j]
              l_int[ind_out, k, i, j] = l_int[ind_out, Nq, i, j]
            end
          end
        end
      end
    end
  end
  nothing
end

function knl_reverse_indefinite_stack_integral!(::Val{dim}, ::Val{N},
                                                ::Val{nvertelem}, auxstate, elems,
                                                ::Val{nout}
                                               ) where {dim, N, nvertelem,
                                                        nout}
  FT = eltype(auxstate)

  Nq = N + 1
  Nqj = dim == 2 ? 1 : Nq

  # note that k is the second not 4th index (since this is scratch memory and k
  # needs to be persistent across threads)
  l_T = MArray{Tuple{nout}, FT}(undef)
  l_V = MArray{Tuple{nout}, FT}(undef)

  @inbounds @loop for eh in (elems; blockIdx().x)
    # Initialize the constant state at zero
    @loop for j in (1:Nqj; threadIdx().y)
      @loop for i in (1:Nq; threadIdx().x)
        ijk = i + Nq * ((j-1) + Nqj * (Nq-1))
        et = nvertelem + (eh - 1) * nvertelem
        @unroll for s = 1:nout
          l_T[s] = auxstate[ijk, s, et]
        end

        # Loop up the stack of elements
        for ev = 1:nvertelem
          e = ev + (eh - 1) * nvertelem
          @unroll for k in 1:Nq
            ijk = i + Nq * ((j-1) + Nqj * (k-1))
            @unroll for s = 1:nout
              l_V[s] = auxstate[ijk, s, e]
            end
            @unroll for s = 1:nout
              auxstate[ijk, nout+s, e] = l_T[s] - l_V[s]
            end
          end
        end
      end
    end
  end
  nothing
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
module ConfidenceIntervals

using Distributions

export binomialCI

# TODO Clopper-Pearson
# TODO Convert to string documentation
# TODO Multinomial CI calculation

## Binomial proportion confidence interval bounds calculation
# Variables
#   - X (int)             = successes
#   - n (int)             = trial
#   - method              = CI calculation method
#   - confidence (float)  = confidence for the interval, 1-α
# Returns
#   lower bound, upper bound
function binomialCI(X::Int, n::Int, method=:wilson, confidence::AbstractFloat=0.95)
  ## Variable checks
  0 <= X || error("Successes must be positive")
  X < n || error("Trials must be larger than successes")
  0 <= confidence <= 1 || error("Confidence must lie between 0 and 1. Typically 0.95")

  ## Method selection
  if method == :agrestiCoull
    return agrestiCoull_CI(X, n, confidence)
  elseif method == :arcsineTransformation
    return arcsineTransformation_CI(X, n, confidence)
  elseif method == :jeffreys
    return jeffreys_CI(X, n, confidence)
  elseif method == :wald
    return wald_CI(X, n, confidence)
  elseif method == :wilson
    return wilson_CI(X, n, confidence)
  elseif method == :wilsonCC
    return wilsonCC_CI(X, n, confidence)
  else
    error("No method found with such name.")
  end
end

function zCalc(confidence)
  return quantile(Normal(), 1-(1-confidence)/2);
end

function agrestiCoull_CI(X, n, confidence)
  ## Variable precalculations
  z = zCalc(confidence)
  m = n+z^2
  μ = (X+z^2/2)/m

  ## Bounds calculation
  σ = sqrt(q*(1-q)/m)
  lower = μ-z*σ
  upper = μ+z*σ

  return lower, upper
end

function arcsineTransformation_CI(X, n, confidence)
  ## Variable precalculations
  p = X/n
  z = zCalc(confidence)

  ## Bounds calculation
  s = z/(2*sqrt(n))
  lower = sin(asin(sqrt(p))-s)^2
  upper = sin(asin(sqrt(p))+s)^2

  return lower, upper
end

function jeffreys_CI(X, n, confidence)
  α = 1-confidence
  lower = quantile(Beta(x+1/2, n-X+1/2), α/2)
  upper = quantile(Beta(x+1/2, n-X+1/2), 1-α/2)
  return lower, upper
end

# Normal approximation
#   - Extremely conservative
function wald_CI(X, n, confidence)
  ## Variable precalculations
  p = X/n
  z = zCalc(confidence)

  ## Bounds calculation
  A = z*sqrt(p*(1-p)/n)
  lower = p-A
  upper = p+A

  return lower, upper
end

# Calculates the Wilson interval bounds, as developed by Edwin Bidwell Wilson (1927, JSTOR 2276774).
#   - Improvement over the normal approximation interval
#   - Good for even small number of trials
#   - Not equal-tailed (systemic bias to centre)
function wilson_CI(X, n, confidence)
  ## Variable precalculations
  neg = n-X
  p = X/n
  z = zCalc(confidence)

  ## Bounds calculation
  #TODO convert to p-X formulas
  A = z*sqrt(X*neg/n+(z^2)/4)
  B = n+z^2
  lower = (X+(z^2)/2-A)/B
  upper = (X+(z^2)/2+A)/B

  return lower, upper
end

# Continuity correction derived from Newcombe (1998, PMID 9595616).
function wilsonCC_CI(X, n, confidence)
  ## Variable precalculations
  p = X/n
  z = zCalc(confidence)

  ## Bounds calculation
  A = 2*n*p+z^2
  B = (z^2)/n+4*n*p*(1-p)
  C = 4*p-2
  D = 2*(n+z^2)
  # Lower bound calculation
  if p == 0
    lower = 0
  else
    lower = maximum([0.0 A-(z*sqrt(B+C)+1)/D])
  end

  # Upper bound calculation
  if p == 1
    upper = 1
  else
    upper = minimum([1.0 A+(z*sqrt(B-C)+1)/D]);
  end

  return lower, upper
end

end

max_stars_repo_name: DCLukas/HypothesisModules

max_issues_repo_name: DCLukas/HypothesisModules

max_forks_repo_name: DCLukas/HypothesisModules


text:
# AUTO GENERATED FILE - DO NOT EDIT

export daq_joystick

"""
    daq_joystick(;kwargs...)

A Joystick component.
A joystick.
Keyword arguments:
- `id` (String; optional): The ID used to identify the color picker in Dash callbacks
- `angle` (Real; optional): Joystick angle in degrees, 0 = right, 90 = up, 180 = left, 270 = down
- `className` (String; optional): Class to apply to the root component element
- `disabled` (Bool; optional): If true, color cannot be picked.
- `force` (Real; optional): Joystick force: distance between cursor and center in big-circle radii
- `label` (optional): Description to be displayed alongside the control. To control styling,
pass an object with label and style properties. label has the following type: String | lists containing elements 'style', 'label'.
Those elements have the following types:
  - `style` (Dict; optional)
  - `label` (String; optional)
- `labelPosition` (a value equal to: 'top', 'bottom'; optional): Where the indicator label is positioned
- `size` (Real; optional): Size (width) of the component in pixels
- `style` (Dict; optional): Style to apply to the root component element
- `theme` (Dict; optional): Theme configuration to be set by a ThemeProvider
"""
function daq_joystick(; kwargs...)
        available_props = Symbol[:id, :angle, :className, :disabled, :force, :label, :labelPosition, :size, :style, :theme]
        wild_props = Symbol[]
        return Component("daq_joystick", "Joystick", "dash_daq", available_props, wild_props; kwargs...)
end