text:
## --------------------------------------------------------------------------
function _extract_sections_to_files_arg_parse(argsv::Vector=ARGS)
    # arg settings
    argset = ArgParse.ArgParseSettings()
    @ArgParse.add_arg_table! argset begin
    "--srcfile", "-f"
        help = "the path to the main file"
        arg_type = String
        default = load_setting("maintex", pwd())
    "--dry-run", "-d"
        help = "Only verbosity, no action"
        action = :store_true
    "--verbose", "-v"
        help = "Print info while executing"
        action = :store_true
    "--overwrite", "-o"
        help = "Allow ovewritting existing files"
        action = :store_true
    "--make-bk", "-b"
        help = "Makes a linear backup"
        action = :store_true
    "--clear-texs", "-c"
        help = "Clear all the tex files excepts the main and the backup"
        action = :store_true
    end
    args = ArgParse.parse_args(argsv, argset)
    to_symbol_arg_dict(args)
end

## --------------------------------------------------------------------------
function extract_sections_to_files(;
        srcfile::String = "",
        overwrite = false,
        dry_run::Bool = false,
        verbose = true,
        make_bk = true,
        names_len = 20,
        clear_texs = true
    )

    # Set up src and dest
    srcdir = abspath(dirname(srcfile))
    !isfile(srcfile) && error(srcfile, " is not a file or is missing!")

    verbose && @info("Args", srcfile, overwrite, dry_run, verbose)
    verbose && println()

    # src lines
    srclines = expand_inputs(;
        srcfile, verbose = false, dry_run = true, commented = false
    )

    # make backup
    if (make_bk && !dry_run)
        bk_lines = map((line) -> string("% ", line), srclines)

        bk_file = joinpath(srcdir, replace(srcfile, ".tex" => ".bk.tex"))
        open(bk_file, "w") do io
            foreach((line) -> println(io, line), bk_lines)
        end
        verbose && @info("Backup created", bk_file)
    end

    # clear
    if (clear_texs && make_bk && !dry_run)
        for file in readdir(srcdir)
            !endswith(file, ".tex") && continue
            endswith(file, ".bk.tex") && continue
            endswith(file, ".expanded.tex") && continue
            file == basename(srcfile) && continue
            path = joinpath(srcdir, file)
            rm(path; force = true)
        end
        verbose && @info("Folder cleared")
    end

    # file tree
    file_tree = Dict{String, Dict}()
    function _register_file!(file, parent)
        haskey(file_tree, file) && error("File name collision: ", file)
        file_tree[file] = Dict(:parent => parent, :lines => String[])
    end

    # main file
    main_file = basename(srcfile)
    _register_file!(main_file, main_file)

    # new files
    sec_counter = 0
    subsec_counter = 0
    subsubsec_counter = 0
    function _newfilename(label)
        fileid = string(sec_counter, ".", subsec_counter, ".", subsubsec_counter, "_")
        label = replace(label, " " => "_")
        if length(label) > names_len
            mid0_end = max(1, div(names_len, 2) - 1)
            mid0 = label[begin:mid0_end]
            mid1_begin = max(mid0_end + 1, div(names_len, 2) + 2)
            mid1 = label[mid1_begin:end]
            label = string(mid0, "...", mid1)
        end
        return string(fileid, label, ".tex")
    end

    # deep record
    MAIN_DEEP = 0
    SECTION_DEEP = 1
    SUBSECTION_DEEP = 2
    SUBSUBSECTION_DEEP = 3
    curr_deep = Dict()
    function _get_parent(deep)
        @assert deep >= 0
        deep -= 1
        while deep >= 0
            haskey(curr_deep, deep) && return curr_deep[deep]
            deep -= 1
        end
        return main_file
    end

    # initi current
    curr_file = main_file
    curr_deep[MAIN_DEEP] = main_file

    # indent
    curr_indent = ""

    for (li, line) in enumerate(srclines)

        # check begin new file
        texsec = extract_tex_section_label(line)
        asssec = extract_texass_sec_name(line)
        texsubsec = extract_tex_subsection_label(line)
        texsubsubsec = extract_tex_subsubsection_label(line)
        isnew_file = !isempty(string(texsec, asssec, texsubsec, texsubsubsec))

        if !isempty(texsec) || !isempty(asssec)
            # update counter
            sec_counter += 1
            subsec_counter = 0
            subsubsec_counter = 0

            new_file = _newfilename(string(texsec, asssec))
            new_deep = SECTION_DEEP

        elseif !isempty(texsubsec)
            # Do subsection stuff
            subsec_counter += 1
            subsubsec_counter = 0

            new_file = _newfilename(texsubsec)
            new_deep = SUBSECTION_DEEP

        elseif !isempty(texsubsubsec)
            # Do subsubsection stuff
            subsubsec_counter += 1

            new_file = _newfilename(texsubsubsec)
            new_deep = SUBSUBSECTION_DEEP

        end
        if isnew_file

            # set new current
            parent = _get_parent(new_deep)
            _register_file!(new_file, parent)
            curr_deep[new_deep] = new_file
            verbose && @info("New file", new_file, parent, li, line)
            curr_file = new_file

            # add input to parent
            curr_indent = extract_indent(line)
            file_dat = file_tree[parent]
            push!(file_dat[:lines],
                string(curr_indent, "\\input{", new_file, "}")
            )

        end

        # main file breakers
        is_begin_document = (extract_tex_begin_label(line) == "document")
        is_begin_document && verbose && @info("Begin document", li, line)
        is_end_document = (extract_tex_end_label(line) == "document")
        is_end_document && verbose && @info("End document", li, line)
        is_end_texass_section = (extract_texass_end_label(line) == "section")
        is_end_texass_section && verbose && @info("TexAss end section", li, line)
        if is_begin_document || is_end_document
            curr_file = main_file
        elseif is_end_texass_section
            # This end goes to the new file
            file_dat = file_tree[curr_file]
            push!(file_dat[:lines], line)
            curr_file = main_file
            continue
        end

        # add line to current file
        file_dat = file_tree[curr_file]
        push!(file_dat[:lines], line)
    end

    ## create files
    function _concat_lines(lines)
        str = join(lines, "\n")
        # small formatting
        str = replace(str, "\n\n\n" => "\n")
        return str
    end

    for (file, content) in file_tree
        path = joinpath(srcdir, file)
        lines = content[:lines]
        if isfile(path) && !overwrite
            verbose && @warn("File already exist", path, overwrite)
        else
            !dry_run && write(path, _concat_lines(lines))
        end
    end

    return file_tree

end


text:
@doc raw"""
    βdifferential_geodesic_startpoint(κ,t,d)

weights for the [`jacobi_field`](@ref) corresponding to the differential of the geodesic
with respect to its start point $D_x g(t;p,q)[X]$.
They are

```math
\beta(\kappa) = \begin{cases}
\frac{\sinh(d(1-t)\sqrt{-\kappa})}{\sinh(d\sqrt{-\kappa})}
&\text{ if }\kappa < 0,\\
1-t & \text{ if } \kappa = 0,\\
\frac{\sin((1-t)d\sqrt{\kappa})}{\sinh(d\sqrt{\kappa})}
&\text{ if }\kappa > 0.
\end{cases}
```

Due to a symmetry agrument, these are also used to compute $D_q g(t; p,q)[\eta]$

# See also

[`differential_geodesic_endpoint`](@ref), [`differential_geodesic_startpoint`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_geodesic_startpoint(κ, t, d)
    (κ == 0) && return 1.0 - t
    (d == 0) && return 1.0 - t
    (κ < 0) && return sinh(sqrt(-κ) * (1.0 - t) * d) / sinh(sqrt(-κ) * d)
    return (κ > 0) && return sin(sqrt(κ) * (1.0 - t) * d) / sin(sqrt(κ) * d)
end
@doc raw"""
    βdifferential_exp_basepoint(κ,t,d)

weights for the [`jacobi_field`](@ref) corresponding to the differential of the geodesic
with respect to its start point $D_p \exp_p X [Y]$. They are

```math
\beta(\kappa) = \begin{cases}
\cosh(\sqrt{-\kappa})&\text{ if }\kappa < 0,\\
1 & \text{ if } \kappa = 0,\\
\cos(\sqrt{\kappa}) &\text{ if }\kappa > 0.
\end{cases}
```

# See also

[`differential_exp_basepoint`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_exp_basepoint(κ, ::Number, d)
    (κ < 0) && return cosh(sqrt(-κ) * d)
    (κ > 0) && return cos(sqrt(κ) * d)
    return 1.0
end
@doc raw"""
    βdifferential_exp_argument(κ,t,d)

weights for the [`jacobi_field`](@ref) corresponding to the differential of the geodesic
with respect to its start point $D_X \exp_p X[Y]$. They are

```math
$\beta(\kappa) = \begin{cases}
\frac{\sinh(d\sqrt{-\kappa})}{d\sqrt{-\kappa}}&\text{ if }\kappa < 0,\\
1 & \text{ if } \kappa = 0,\\
\frac{\sin(d\sqrt{\kappa})}{\sqrt{d\kappa}}&\text{ if }\kappa > 0.
\end{cases}
```

# See also

[`differential_exp_argument`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_exp_argument(κ, ::Number, d)
    ((κ == 0) || (d == 0)) && return 1.0
    (κ < 0) && return sinh(sqrt(-κ) * d) / (d * sqrt((-κ)))
    return (κ > 0) && return sin(sqrt(κ) * d) / (d * sqrt(κ))
end
@doc raw"""
    βdifferential_log_basepoint(κ,t,d)

weights for the [`jacobi_field`](@ref) corresponding to the differential of the geodesic
with respect to its start point $D_p \log_p q[X]$. They are

```math
\beta(\kappa) = \begin{cases}
-\sqrt{-\kappa}d\frac{\cosh(d\sqrt{-\kappa})}{\sinh(d\sqrt{-\kappa})}&\text{ if }\kappa < 0,\\
-1 & \text{ if } \kappa = 0,\\
-\sqrt{\kappa}d\frac{\cos(d\sqrt{\kappa})}{\sin(d\sqrt{\kappa})}&\text{ if }\kappa > 0.
\end{cases}
```

# See also

[`differential_log_argument`](@ref), [`differential_log_argument`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_log_basepoint(κ, ::Number, d)
    ((d == 0) || (κ == 0)) && return -1.0
    (κ < 0) && return -sqrt(-κ) * d * cosh(sqrt(-κ) * d) / sinh(sqrt(-κ) * d)
    return (κ > 0) && return -sqrt(κ) * d * cos(sqrt(κ) * d) / sin(sqrt(κ) * d)
end
@doc raw"""
    βdifferential_log_argument(κ,t,d)

weights for the JacobiField corresponding to the differential of the logarithmic
map with respect to its argument $D_q \log_p q[X]$. They are

```math
\beta(\kappa) = \begin{cases}
\frac{ d\sqrt{-\kappa} }{\sinh(d\sqrt{-\kappa})}&\text{ if }\kappa < 0,\\
1 & \text{ if } \kappa = 0,\\
\frac{ d\sqrt{\kappa} }{\sin(d\sqrt{\kappa})}&\text{ if }\kappa > 0.
\end{cases}
```

# See also

[`differential_log_basepoint`](@ref), [`jacobi_field`](@ref)
"""
function βdifferential_log_argument(κ, ::Number, d)
    ((d == 0) || (κ == 0)) && return 1.0
    (κ < 0) && return sqrt(-κ) * d / sinh(sqrt(-κ) * d)
    return sqrt(κ) * d / sin(sqrt(κ) * d)
end

@doc raw"""
    Y = adjoint_Jacobi_field(M, p, q, t, X, β)

Compute the AdjointJacobiField $J$ along the geodesic $γ_{p,q}$ on the manifold
$\mathcal M$ with initial conditions (depending on the application)
$X ∈ T_{γ_{p,q}(t)}\mathcal M$ and weights $β$. The result is a vector
$Y ∈ T_p\mathcal M$. The main difference to [`jacobi_field`](@ref) is the,
that the input `X` and the output `Y` switched tangent spaces.
For detais see [`jacobi_field`](@ref)
"""
function adjoint_Jacobi_field(M::Manifold, p, q, t, X, β=βdifferential_geodesic_startpoint)
    x = shortest_geodesic(M, p, q, t)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M, p, q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.(Ref(M), Ref(p), V, Ref(x), Ref(ParallelTransport()))
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y = sum(
        (inner.(Ref(M), Ref(x), Ref(X), Θ)) .*
        (β.(B.data.eigenvalues, Ref(t), Ref(distance(M, p, q)))) .* V,
    )
    return Y
end
function adjoint_Jacobi_field(
    M::Circle{ℝ}, p::Real, q::Real, t::Real, X::Real, β=βdifferential_geodesic_startpoint
)
    x = shortest_geodesic(M, p, q, t)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M, p, q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.(Ref(M), Ref(p), V, Ref(x), Ref(ParallelTransport()))
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y = sum(
        (inner.(Ref(M), Ref(x), Ref(X), Θ)) .*
        (β.(B.data.eigenvalues, Ref(t), Ref(distance(M, p, q)))) .* V,
    )[1]
    return Y
end
function adjoint_Jacobi_field(
    M::AbstractPowerManifold, p, q, t, X, β=βdifferential_geodesic_startpoint
)
    rep_size = representation_size(M.manifold)
    Y = allocate_result(M, adjoint_Jacobi_field, p, X)
    for i in get_iterator(M)
        #lY = adjoint_Jacobi_field(M.manifold, p[i], q[i], t, X[i], β)
        Y[M, i] = adjoint_Jacobi_field(M.manifold, p[M, i], q[M, i], t, X[M, i], β)
    end
    return Y
end

@doc raw"""
    Y = jacobi_field(M, p, q, t, X, β)
compute the Jacobi jield $J$ along the geodesic $γ_{p,q}$ on the manifold $\mathcal M$ with
initial conditions (depending on the application) $X ∈ T_p\mathcal M$ and weights $β$. The
result is a tangent vector `Y` from $T_{γ_{p,q}(t)}\mathcal M$.

# See also

[`adjoint_Jacobi_field`](@ref)
"""
function jacobi_field(M::Manifold, p, q, t, X, β=βdifferential_geodesic_startpoint)
    x = shortest_geodesic(M, p, q, t)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M, p, q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.(Ref(M), Ref(p), V, Ref(x), Ref(ParallelTransport()))
    Y = zero_tangent_vector(M, p)
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y .= sum(
        (inner.(Ref(M), Ref(p), Ref(X), V)) .*
        (β.(B.data.eigenvalues, Ref(t), Ref(distance(M, p, q)))) .* Θ,
    )
    return Y
end
function jacobi_field(
    M::Circle{ℝ}, p::Real, q::Real, t::Real, X::Real, β=βdifferential_geodesic_startpoint
)
    x = shortest_geodesic(M, p, q, t)
    B = get_basis(M, p, DiagonalizingOrthonormalBasis(log(M, p, q)))
    V = get_vectors(M, p, B)
    Θ = vector_transport_to.(Ref(M), Ref(p), V, Ref(x), Ref(ParallelTransport()))
    Y = zero_tangent_vector(M, p)
    # Decompose wrt. frame, multiply with the weights from w and recompose with Θ.
    Y = sum(
        (inner.(Ref(M), Ref(p), Ref(X), V)) .*
        (β.(B.data.eigenvalues, Ref(t), Ref(distance(M, p, q)))) .* Θ,
    )[1]
    return Y
end
function jacobi_field(
    M::AbstractPowerManifold, p, q, t, X, β=βdifferential_geodesic_startpoint
)
    rep_size = representation_size(M.manifold)
    Y = allocate_result(M, adjoint_Jacobi_field, p, X)
    for i in get_iterator(M)
        Y[M, i] = jacobi_field(M.manifold, p[M, i], q[M, i], t, X[M, i], β)
    end
    return Y
end

ext: jl


text:
using Cassette
using SpecializeVarargs: @specialize_vararg
using ChainRules
using ChainRulesCore
import ChainRulesCore: Zero

using Cassette: overdub, Context, nametype, similarcontext

Cassette.@context DualContext

const TaggedCtx{T} = Context{nametype(DualContext),T}

function dualcontext()
    # Note that the `dualtag()` is not of the same type as that of the
    # Duals constructed in this context, because it is called in the older context
    # New duals are of type Tag{typeof(dualtag())}
    Cassette.disablehooks(DualContext(metadata=dualtag(), pass=CustomDispatchPass))
end

# Calls to `dualtag` are aware of the current context. Note that the tags
# produced in the current context are `Tag{T} where T` is the metadata type of
# the context.
@inline Cassette.overdub(ctx::TaggedCtx{T}, ::typeof(dualtag)) where {T} = Tag{T}()

@inline Cassette.overdub(ctx::TaggedCtx, ::typeof(find_dual), args...) = find_dual(args...)

@inline _value(::Any, x) = x
@inline _value(::Tag{T}, d::Dual{Tag{T}}) where T = value(d)

@inline _partials(::Any, x) = Zero()
@inline _partials(::Tag{T}, d::Dual{Tag{T}}) where T = d.partials

@inline _values(S, xs) = map(x->_value(S, x), xs)
@inline _partialss(S, xs) = map(x->_partials(S, x), xs)

@inline anydual(x, xs...) = anydual(xs...)
@inline anydual(x::Dual, xs...) = true
@inline anydual() = false

# necessary special cases:

# this makes splatting work...
@inline isinteresting(ctx::TaggedCtx, f::typeof(Core._apply), g, xs...) = Core._apply(isinteresting, (ctx, g), xs...)
@inline function alternative(ctx::TaggedCtx{T}, f::typeof(Core._apply), g, args...) where {T}
    Core._apply(alternative, (ctx, g), args...)
end

# actually interesting:
@inline isinteresting(ctx::TaggedCtx, f, a) = anydual(a)
@inline isinteresting(ctx::TaggedCtx, f, a, b) = anydual(a, b)
@inline isinteresting(ctx::TaggedCtx, f, a, b, c) = anydual(a, b, c)
@inline isinteresting(ctx::TaggedCtx, f, a, b, c, d) = anydual(a, b, c, d)
@inline isinteresting(ctx::TaggedCtx, f, args...) = anydual(args...)
@inline isinteresting(ctx::TaggedCtx, f::Core.Builtin, args...) = false
@inline isinteresting(ctx::TaggedCtx, f::Union{typeof(ForwardDiff2.find_dual),
                                               typeof(ForwardDiff2.anydual)}, args...) = false

@specialize_vararg 4 @inline function _frule_overdub2(ctx::TaggedCtx{T}, f::F, args...) where {T,F}
    # Here we can assume that one or more `args` is a Dual with tag
    # of type T.

    tag = Tag{T}()
    # unwrap only duals with the tag T.
    vs = _values(tag, args)

    # extract the partials only for the current tag
    # so we can pass them to the pushforward
    ps = _partialss(tag, args)

    # default `dself` to `Zero()`
    dself = Zero()

    # call frule to see if there is a rule for this call:
    if ctx.metadata isa Tag
        ctx1 = similarcontext(ctx, metadata=oldertag(ctx.metadata))

        # we call frule with an older context because the Dual numbers may
        # themselves contain Dual numbers that were created in an older context
        frule_result = overdub(ctx1, frule, f, vs..., dself, ps...)
    else
        frule_result = frule(f, vs..., dself, ps...)
    end

    if frule_result === nothing
        # this means there is no frule
        # We can't just do f(args...) here because `f` might be
        # a closure which closes over a Dual number, hence we call
        # recurse. Recurse overdubs the calls inside `f` and not `f` itself
        return Cassette.overdub(ctx, f, args...)
    else
        # this means there exists an frule for this specific call.
        # frule_result is then a tuple (val, pushforward) where val
        # is the primal result. (Note: this may be Dual numbers but only
        # with an older tag)
        val, ∂s = frule_result

        return if ∂s isa Tuple
            map(val, ∂s) do v, ∂
                ∂ isa Zero ? v : Dual{Tag{T}}(v, ∂)
            end
        else
            ∂s isa Zero ? val : Dual{Tag{T}}(val, ∂s)
        end
    end
end

@specialize_vararg 4 @inline function alternative(ctx::TaggedCtx{T}, f::F, args...) where {T,F}
    # This method only executes if `args` contains at least 1 Dual
    # the question is what is its tag

    tag = Tag{T}()

    idx = find_dual(tag, args...)
    if f === Dual
        return overdub(ctx, f, args...)
    elseif idx === 0
        # This is the base case for the recursion in this function which
        # tries to do the alternative with successively older contexts
        # the oldest context is not a `Tag`
        if !(T<:Tag)
            return Cassette.fallback(ctx, f, args...)
        end
        # none of the arguments have the same tag as the context
        # try with the parent context
        ctx1 = similarcontext(ctx, metadata=oldertag(ctx.metadata))
        return overdub(ctx1, ()->f(args...)) # This will wrap it in isinteresting; alternative
    else
        # call ChainRules.frule to execute `f` and
        # get a function that computes the partials
        return _frule_overdub2(ctx, f, args...)
    end
end

function dualrun(f, args...)
    ctx = dualcontext()
    return overdub(ctx, f, args...)
end

const BINARY_PREDICATES = Symbol[:isequal, :isless, :<, :>, :(==), :(!=), :(<=), :(>=)]

for pred in BINARY_PREDICATES
    @eval function alternative(ctx::TaggedCtx, ::typeof($(pred)), x, y)
        vx, vy = value(x), value(y)
        return isinteresting(ctx, $pred, vx, vy) ?
        alternative(ctx, $pred, vx, vy) :
        $pred(vx, vy)
    end
end

#### recursion early termination condition
@inline Cassette.overdub(ctx::TaggedCtx, f::Core.Builtin, args...) = f(args...)
@inline Cassette.overdub(ctx::TaggedCtx{T}, f::Union{typeof(value),typeof(partials)}, d::Dual{T}) where {T<:Tag{Nothing}} = f(d)
@inline Cassette.overdub(ctx::TaggedCtx{T}, f::typeof(allpartials), d::AbstractArray{<:Dual{T}})  where {T<:Tag{Nothing}} = f(d)
@inline Cassette.overdub(ctx::TaggedCtx, f::typeof(Base.Broadcast.check_broadcast_axes), args...) = f(args...)
@inline Cassette.overdub(ctx::TaggedCtx, f::typeof(Base.Broadcast.preprocess_args), args...) = f(args...)

##### Inference Hacks
@inline isinteresting(ctx::TaggedCtx, f::Union{typeof(Base.print_to_string),typeof(hash)}, args...) = false
@noinline Cassette.overdub(ctx::TaggedCtx, f::Union{typeof(Base.print_to_string),typeof(hash),typeof(Core.throw)}, args...) = f(args...)


text:
using CUDA, Test, AutoGrad, Knet
using CUDA.CUDNN:
    cudnnActivationForward,
    cudnnActivationForward!,
    cudnnActivationBackward,
    cudnnActivationDescriptor,
        cudnnActivationDescriptor_t,
        cudnnCreateActivationDescriptor,
        cudnnSetActivationDescriptor,
        cudnnGetActivationDescriptor,
        cudnnDestroyActivationDescriptor,
    cudnnActivationMode_t,
        CUDNN_ACTIVATION_SIGMOID,
        CUDNN_ACTIVATION_RELU,
        CUDNN_ACTIVATION_TANH,
        CUDNN_ACTIVATION_CLIPPED_RELU,
        CUDNN_ACTIVATION_ELU,
        CUDNN_ACTIVATION_IDENTITY,
    cudnnNanPropagation_t,
        CUDNN_NOT_PROPAGATE_NAN,
        CUDNN_PROPAGATE_NAN,
    handle

@testset "cudnn/activation" begin

    x = Param(CUDA.randn(Float64,10))
    @test @gcheck cudnnActivationForward(x, alpha=1, coef=1, nanOpt=CUDNN_NOT_PROPAGATE_NAN, mode=CUDNN_ACTIVATION_SIGMOID)
    @test @gcheck cudnnActivationForward(x, alpha=1, coef=1, nanOpt=CUDNN_NOT_PROPAGATE_NAN, mode=CUDNN_ACTIVATION_RELU)
    @test @gcheck cudnnActivationForward(x, alpha=1, coef=1, nanOpt=CUDNN_NOT_PROPAGATE_NAN, mode=CUDNN_ACTIVATION_TANH)
    @test @gcheck cudnnActivationForward(x, alpha=1, coef=1, nanOpt=CUDNN_NOT_PROPAGATE_NAN, mode=CUDNN_ACTIVATION_CLIPPED_RELU)
    @test @gcheck cudnnActivationForward(x, alpha=1, coef=1, nanOpt=CUDNN_NOT_PROPAGATE_NAN, mode=CUDNN_ACTIVATION_ELU)
    @test @gcheck cudnnActivationForward(x, alpha=2, coef=1, nanOpt=CUDNN_NOT_PROPAGATE_NAN, mode=CUDNN_ACTIVATION_SIGMOID)
    @test @gcheck cudnnActivationForward(x, alpha=1, coef=2, nanOpt=CUDNN_NOT_PROPAGATE_NAN, mode=CUDNN_ACTIVATION_CLIPPED_RELU)
    @test @gcheck cudnnActivationForward(x, alpha=1, coef=2, nanOpt=CUDNN_NOT_PROPAGATE_NAN, mode=CUDNN_ACTIVATION_ELU)
    @test @gcheck cudnnActivationForward(x, alpha=1, coef=1, nanOpt=CUDNN_PROPAGATE_NAN,     mode=CUDNN_ACTIVATION_SIGMOID)

end


text:
using StructArrays # for type definitions
using Statistics # for random input generation
using BenchmarkTools # for benchmark
# using DataFrames
using GFlops # to count gflops
using Suppressor

export Funb, FunbArray, benchmark!

################################################################
"""
    Funb(;fun, limits, types, dims)

Creates random inputs for a function based on limits, types, and dims specified.

# Arguments
- fun: the function `:fun` or :(Module.fun)
- limits: min and max of possible values
- types : type of elements
- dims: Array of dimensions of the input vectors for each argument. Each column is for a new set of sizes, and each row is for different input arguments. So:
     - each element gives the size of the input, and it is a:
        - Number (for 1D)
        - Tuple (for N-D)
     - each row for each function argument
     - each column for each dimension set

# Examples
```julia
config = Funb(
    fun = sin,
    limits = [(-1,1)],
    types = [Float32, Float64],
    dims = [100],
)
```
or just in a compact form:
```julia
config = Funb( sin, [(-1,1)], [Float32, Float64], [100])
```
"""
struct Funb

    # public
    fun::T where {T<:Function}
    limits::Vector{Tuple{T1,T2}} where {T1,T2}
    types::Vector{DataType}
    dims::VecOrMat{T1} where {T1<:Union{T2, Tuple}} where {T2<:Number}

    # private - calculated

    sets #::Vector{Vector{T}} where {T<:Distribution} # vector of input sets for each function. Each set is an array of distributions
    inputs #::Vector{Vector{T}} where {T}
    results
    median

    gflops
end

Funb( fun, limits, types, dims) =  Funb( fun = fun, limits = limits, types = types, dims = dims)

function Funb(; fun, limits, types, dims)

        numTypes  = length(types)
        sets = Vector(undef, numTypes)
        inputs = Vector(undef, numTypes)
        results = Vector{Any}(undef, numTypes)
        median = Vector{Any}(undef, numTypes)
        gflops = Vector{Any}(undef, numTypes)

        for iType = 1:numTypes

            numArgs, numDimsSets = numArgsDims(dims)

            sets[iType] = Vector(undef, numArgs)
            inputs[iType] = Vector(undef, numDimsSets)

            results[iType] = Vector{BenchmarkTools.Trial}(undef, numDimsSets)
            median[iType] = Vector{Float64}(undef, numDimsSets)
            gflops[iType] = Vector{Float64}(undef, numDimsSets)

            for iDimSet = 1:numDimsSets

                inputs[iType][iDimSet] = Vector(undef, numArgs) # {Array{types[iType]}}

                for iArg = 1:numArgs
                    try
                        sets[iType][iArg] = Uniform2(types[iType], limits[iArg]...)
                        inputs[iType][iDimSet][iArg] = Base.rand(sets[iType][iArg], dims[iArg, iDimSet])
                    catch
                        sets[iType][iArg] = Uniform(types[iType], limits[iArg]...)
                        inputs[iType][iDimSet][iArg] = Base.rand(sets[iType][iArg], dims[iArg, iDimSet])
                    end
                end
            end

        end

    return Funb( fun, limits, types, dims, sets, inputs, results, median, gflops)
end
################################################################
"""
    FunbArray

Array of Funb configs for different functions.

# Examples
```julia
using AcuteBenchmark

configs = FunbArray([
    Funb( sin, [(-1,1)],[Float32, Float64], [100] );
    Funb( atan, [(-1,1), (-1,1)],[Float32, Float64],[100, 100] );
    Funb( *, [(-1, 1), (-1, 1), (-1, 1)], [Float32, Float64], [(100,100), (100,100)] );
    ])

```

You can also directly give the configs in vectors:
```julia
configs = FunbArray(
    fun =   [sin,
             atan,
             *],
    limits = [[(-1,1)],
             [(-1,1), (-1,1)],
             [(-1, 1), (-1, 1), (-1, 1)]],
    types = fill([Float32, Float64], (3)),
    dims = [ [100],
             [100, 100],
             [(100,100), (100,100)] ],
)
```
"""
FunbArray(configs::Array{Funb}) = StructArray(configs)
FunbArray(configs::Funb) = StructArray([configs])

function FunbArray(;fun, limits, types, dims)

    numFun = length(fun)
    configs = Vector{Funb}(undef, numFun)

    for iFun = 1:numFun
        configs[iFun] = Funb(fun = fun[iFun], limits = limits[iFun], types = types[iFun], dims = dims[iFun])
    end

    return StructArray(configs)
end
################################################################
"""
    benchmark!(config::StructArray{Funb}) # FunbArray{Funb}
    benchmark!(config::Array{Funb})

Performs the benchmarking on a given Funb.

# Examples
```julia
using AcuteBenchmark

configs = FunbArray([
    Funb( sin, [(-1,1)],[Float32, Float64], [100] );
    Funb( atan, [(-1,1), (-1,1)],[Float32, Float64],[100, 100] );
    Funb( *, [(-1, 1), (-1, 1), (-1, 1)], [Float32, Float64], [(100,100), (100,100)] );
    ])

benchmark!(configs)
```
"""
function benchmark!(config::StructArray{Funb})

    numFuns = length(config.fun)

    for (iFun, fun) in enumerate(config.fun)

        numTypes = length(config.types[iFun])

        for (iType, type) in enumerate(config.types[iFun])

            numArgs, numDimsSets = numArgsDims(config.dims[iFun])

            for iDimSet = 1:numDimsSets
                inp = config.inputs[iFun][iType][iDimSet]

                println("Benchmarking $fun - $type - dimension set $iDimSet")

                if numArgs == 1 # single argument function
                    if hasmethod(fun, (typeof(inp[1]),)) # check if array method exists
                        config[iFun].results[iType][iDimSet] = @benchmark $fun($inp[1])

                        # counting gflops
                        @suppress config[iFun].gflops[iType][iDimSet] = @gflops $fun($inp[1])

                    else # broadcast
                        config[iFun].results[iType][iDimSet] = @benchmark $fun.($inp[1])

                        # counting gflops - TODO
                        @suppress config[iFun].gflops[iType][iDimSet] = length(inp[1]) * @gflops $fun($inp[1][1])
                    end
                else
                    if hasmethod(fun, Tuple(typeof.(inp))) # check if array method exists
                        config[iFun].results[iType][iDimSet] = @benchmark $fun($inp...)

                        # counting gflops
                        @suppress config[iFun].gflops[iType][iDimSet] = @gflops $fun($inp...)
                    else # broadcast
                        config[iFun].results[iType][iDimSet] = @benchmark $fun.($inp...)

                        # counting gflops - TODO
                        inp1dim = zeros(config[iFun].types[iType], numArgs)
                        for iArg = 1:numArgs
                            inp1dim[iArg] = Base.rand(config[iFun].sets[iType][iArg])
                        end
                        @suppress config[iFun].gflops[iType][iDimSet] = length(inp[1]) * @gflops $fun($inp1dim...)
                    end
                end
                config[iFun].median[iType][iDimSet] = median(config[iFun].results[iType][iDimSet].times) / 1000 # micro seconds

            end
        end
    end
    return config
end

benchmark!(config::Union{Funb, Array{Funb}}) = benchmark!(FunbArray(config))


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    components(adjacency)

Return connected components of `adjacency`.
"""
function components(adjacency::AbstractMatrix{Int})
  n = size(adjacency, 1)
  vertices = Set(1:n)
  components = Vector{Int}[]
  while !isempty(vertices)
    u = pop!(vertices)
    c = component(adjacency, u)
    for v in setdiff(c, [u])
      pop!(vertices, v)
    end
    push!(components, c)
  end
  components
end

"""
    component(adjacency, vertex)

Return connected component of `adjacency` containing `vertex.
"""
function component(adjacency::AbstractMatrix{Int}, vertex::Int)
  frontier = [vertex]
  visited  = Int[]
  # breadth-first search
  while !isempty(frontier)
    u = pop!(frontier)
    push!(visited, u)
    for v in findall(!iszero, adjacency[u,:])
      if v ∉ visited
        push!(frontier, v)
      end
    end
  end
  visited
end

"""
    isacyclic(adjacency)

Check whether or not the `adjacency` matrix is acyclic
using Kahn's topological sort.
"""
function isacyclic(adjacency::AbstractMatrix{Int})
  # copy input to avoid side effects
  A = copy(adjacency)

  # find root nodes
  roots = Vector{Int}()
  for j in 1:size(A, 2)
    if all(A[:,j] .== 0)
      push!(roots, j)
    end
  end

  # Kahn's algorithm
  sorted = Vector{Int}()
  while !isempty(roots)
    i = pop!(roots)
    push!(sorted, i)

    # for all edges i → j
    for j in findall(A[i,:] .== 1)
      # remove edge i → j
      A[i,j] = 0

      # if j has no other incoming edge
      if all(A[:,j] .== 0)
        push!(sorted, j)
      end
    end
  end

  any(A .== 1) ? false : true
end


text:
# This routine ABSOLUTELY should not allocate any memory, called inside time loop.
@inbounds @fastmath function compute_gradient!(issp::Int64, pac::PFdtdc, pass::Vector{PFdtdss}, pap::PFdtdp)
	# aliases
	p=pap.p
	pp=pap.pp
	ppp=pap.ppp
	dpdx=pap.dpdx
	dpdz=pap.dpdz
	δtI=pac.δtI
	grad_modtt=pass[issp].grad_modtt
	grad_modrrvx=pass[issp].grad_modrrvx
	grad_modrrvz=pass[issp].grad_modrrvz

	gmodtt!(grad_modtt,p,pp,ppp,pac.nx,pac.nz,δtI)
	gmodrrvx!(grad_modrrvx,dpdx,pac.nx,pac.nz)
	gmodrrvz!(grad_modrrvz,dpdz,pac.nx,pac.nz)
end

@inbounds @fastmath function gmodtt!(grad_modtt,p,pp,ppp,nx,nz,δtI,)
	pppw=ppp[1]
	pw=p[1]
	ppw=pp[1]
	ppw2=pp[2]
	for ix=1:nx
		@simd for iz=1:nz
			# p at [it], pp at [it-1]	# dpdx and dpdz at [it]	# gradients w.r.t inverse of modtt, i.e., 1/rho/c^2
			@inbounds grad_modtt[iz,ix] += ((-1.0 * (pppw[iz, ix, 1] + pw[iz, ix,  1] - 2.0 * ppw[iz, ix,  1]) * δtI * δtI) *  ppw2[iz,ix,1])
		end
	end
end

@inbounds @fastmath function gmodrrvx!(grad_modrrvx,dpdx,nx,nz)
	dpdxw1=dpdx[1]
	dpdxw2=dpdx[2]
	for ix=1:nx
		@simd for iz=1:nz
			# gradient w.r.t. inverse of rho on vx and vz grids
			@inbounds grad_modrrvx[iz,ix] += (- dpdxw2[iz,ix,1]*dpdxw1[iz,ix,1])
		end
	end

end
@inbounds @fastmath function gmodrrvz!(grad_modrrvz,dpdz,nx,nz)
	dpdzw1=dpdz[1]
	dpdzw2=dpdz[2]
	for ix=1:nx
		@simd for iz=1:nz
			# gradient w.r.t. inverse of rho on vx and vz grids
			@inbounds grad_modrrvz[iz,ix] += (- dpdzw2[iz,ix,1]*dpdzw1[iz,ix,1])
		end
	end
end

@inbounds @fastmath function scale_gradient!(issp::Int64,pass::Vector{PFdtdss},δ)
	grad_modtt=pass[issp].grad_modtt
	grad_modrrvx=pass[issp].grad_modrrvx
	grad_modrrvz=pass[issp].grad_modrrvz
	"gradient is formed by intergration over time, hence multiply with δt, but why not?"
	"I don't completely understand where the factors δx and δz are coming from..."
	"probably the source term should not be multiplied by δxI and δzI during adjoint propagation"
	rmul!(grad_modtt,δ)
	rmul!(grad_modrrvx,δ)
	rmul!(grad_modrrvz,δ)
end

function grad_modrr!(pac::PFdtdc)
	grad_modrr_sprayrr!(pac.grad_modrr_stack,pac.grad_modrrvx_stack,pac.grad_modrrvz_stack)
	grad_modrr_sprayrrvx!(pac.grad_modrr_stack,pac.grad_modrrvx_stack)
	grad_modrr_sprayrrvz!(pac.grad_modrr_stack,pac.grad_modrrvz_stack)
end

function stack_grads!(pac::PFdtdc, pap::PFdtdp)
	# theses are SharedArrays
	gmodtt=pac.grad_modtt_stack
	gmodrrvx=pac.grad_modrrvx_stack
	gmodrrvz=pac.grad_modrrvz_stack
	pass=pap.ss
	for issp in 1:length(pass)
		gs=pass[issp].grad_modtt
		for i in eachindex(gmodtt)
			gmodtt[i]+=gs[i]  # only update gmodtt
		end
		gs=pass[issp].grad_modrrvx
		for i in eachindex(gmodrrvx)
			gmodrrvx[i]+=gs[i]
		end
		gs=pass[issp].grad_modrrvz
		for i in eachindex(gmodrrvz)
			gmodrrvz[i]+=gs[i]
		end
	end
end

function update_gradient!(pac::PFdtdc)
	nx, nz=pac.nx, pac.nz
	nznxd = prod(length.(pac.model.mgrid))

	# combine rrvx and rrvz
	grad_modrr!(pac)

	gradient=pac.gradient
	# truncate
	gmodtt=view(pac.grad_modtt_stack,npml+1:nz-npml,npml+1:nx-npml)
	gmodrr=view(pac.grad_modrr_stack,npml+1:nz-npml,npml+1:nx-npml)
	for i in 1:nznxd
		# parameterization is  [:KI, :ρI, :null]
		gradient[i]=gmodtt[i]  # update gmodtt
		gradient[nznxd+i]=gmodrr[i]  # update gmodrr
	end
end


text:
using IterativeSolvers
using FactCheck

facts("Basic operations") do

context("Adivtype") do
    A = rand(1:10, 5, 5)
    b = rand(Float32, 5)
    @fact IterativeSolvers.Adivtype(A, b) --> Float32
end

context("Linear operator defined as a function") do
    # A = cycle-back operator
    function shiftback!(output, b)
        n = length(b)
        length(output) == n || throw(DimensionMismatch())
        for i = 1:n-1
            output[i+1] = b[i]
        end
        output[1] = b[n]
        output
    end

    # A' = cycle-forward operator
    function shiftfwd!(output, b)
        n = length(b)
        length(output) == n || throw(DimensionMismatch())
        for i = 2:n
            output[i-1] = b[i]
        end
        output[n] = b[1]
        output
    end

    b = rand(Float32, 5)
    Atimesb = [b[end]; b[1:end-1]]
    Aptimesb = [b[2:end]; b[1]]

    A = MatrixFcn{Int}(5, 5, shiftback!)
    @fact eltype(A) --> Int
    @fact size(A) --> (5,5)
    @fact size(A,1) --> 5
    @fact size(A,2) --> 5
    @fact length(A) --> 25
    @fact A*b --> Atimesb

    output = similar(b)
    @fact A_mul_B!(output, A, b) --> Atimesb
    @fact A'*b --> Atimesb

    A = MatrixCFcn{Int}(5, 5, shiftback!, shiftfwd!)
    @fact eltype(A) --> Int
    @fact size(A) --> (5,5)
    @fact size(A,1) --> 5
    @fact size(A,2) --> 5
    @fact length(A) --> 25
    @fact A*b --> Atimesb

    output = similar(b)
    @fact A_mul_B!(output, A, b) --> Atimesb
    @fact A'*b --> Aptimesb
    @fact Ac_mul_B!(output, A, b) --> Aptimesb
end

end


text:
# Decompressor Codec
# ==================

struct XzDecompressor <: TranscodingStreams.Codec
    stream::LZMAStream
    memlimit::Integer
    flags::UInt32
end

function Base.show(io::IO, codec::XzDecompressor)
    print(io, summary(codec), "(memlimit=$(codec.memlimit), flags=$(codec.flags))")
end

const DEFAULT_MEM_LIMIT = typemax(UInt64)

"""
    XzDecompressor(;memlimit=$(DEFAULT_MEM_LIMIT), flags=LZMA_CONCATENATED)

Create an xz decompression codec.

Arguments
---------
- `memlimit`: memory usage limit as bytes
- `flags`: decoder flags
"""
function XzDecompressor(;memlimit::Integer=DEFAULT_MEM_LIMIT, flags::UInt32=LZMA_CONCATENATED)
    if memlimit ≤ 0
        throw(ArgumentError("memlimit must be positive"))
    end
    # NOTE: flags are checked in liblzma
    return XzDecompressor(LZMAStream(), memlimit, flags)
end

const XzDecompressorStream{S} = TranscodingStream{XzDecompressor,S} where S<:IO

"""
    XzDecompressorStream(stream::IO; kwargs...)

Create an xz decompression stream (see `XzDecompressor` for `kwargs`).
"""
function XzDecompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:memlimit, :flags))
    return TranscodingStream(XzDecompressor(;x...), stream; y...)
end

# Methods
# -------

function TranscodingStreams.initialize(codec::XzDecompressor)
    ret = stream_decoder(codec.stream, codec.memlimit, codec.flags)
    if ret != LZMA_OK
        lzmaerror(codec.stream, ret)
    end
    return
end

function TranscodingStreams.finalize(codec::XzDecompressor)
    free(codec.stream)
end

function TranscodingStreams.startproc(codec::XzDecompressor, mode::Symbol, error::Error)
    ret = stream_decoder(codec.stream, codec.memlimit, codec.flags)
    if ret != LZMA_OK
        error[] = ErrorException("xz error")
        return :error
    end
    return :ok
end

function TranscodingStreams.process(codec::XzDecompressor, input::Memory, output::Memory, error::Error)
    stream = codec.stream
    stream.next_in = input.ptr
    stream.avail_in = input.size
    stream.next_out = output.ptr
    stream.avail_out = output.size
    if codec.flags & LZMA_CONCATENATED != 0
        action = stream.avail_in > 0 ? LZMA_RUN : LZMA_FINISH
    else
        action = LZMA_RUN
    end
    ret = code(stream, action)
    Δin = Int(input.size - stream.avail_in)
    Δout = Int(output.size - stream.avail_out)
    if ret == LZMA_OK
        return Δin, Δout, :ok
    elseif ret == LZMA_STREAM_END
        return Δin, Δout, :end
    else
        error[] = ErrorException(lzma_error_string(ret))
        return Δin, Δout, :error
    end
end


text:
using Profile
using Dorothy
using Formats

profilepath = @__DIR__
datapath = joinpath(@__DIR__, "..", "..", "data")

readf("$(datapath)/PLC.pdb")

@profile (for i = 1:2; readf("$(datapath)/PLC.pdb"); end)

open("$(profilepath)/pdbread.txt", "w") do io
	Profile.print(IOContext(io, :displaysize => (24, 500)))
end


text:
using MemorableUniqueIdentifier, Test

@testset "basic" begin
    @test create(6) !== nothing
    @test animal("1200000000000000000000010c6228c1") === "Gloomless Fly"
    @test search("fe3a1eb0bca7542150e37ce4022a366b") === "Female Bobcat"
    @test validate("1200000000000000000000010c6228c1") === true
    @test difficulty("1200000000000000000000010c6228c1") === 12
    @test length(mine_until(6, 3)) === 3
end

max_stars_repo_name: rustyconover/muid-julia

max_issues_repo_name: rustyconover/muid-julia

max_forks_repo_name: rustyconover/MemorableUniqueIdentifier


text:
module SeenVars

using ..Traits, ..Tagging, ..Var
import ..Var

export tagseen

tagseen(ω, seen = Set()) = tag(ω, (seen = seen,))

function Var.posthook(::trait(Seen), ret, f::PrimRandVar, ω)
  println("does this handle i.i.d right?")
  push!(ω.tags.seen, f)
end

export ctxapply
function ctxapply(x, ω)
  ..
end

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# Big TODO we are still missing the functionality that does the aligning in the constructors

Assemblage(assm::SEAssemblage) = Assemblage(assm.site, assm.occ) # Not a copy constructor - Just a function that will reduce derived types to the base type

# a constructor that takes occ and coords as one single DataFrame format and separates them
function Assemblage(occ::DataFrames.DataFrame; kwargs...)
  occ, coords = parsesingleDataFrame(occ)
  Assemblage(occ, coords; kwargs...)
end

# a constructor that takes occ as a DataFrame
Assemblage(occ::DataFrames.DataFrame, coords::Union{AbstractMatrix, DataFrames.DataFrame}; sitecolumns = true, kwargs...) = Assemblage(ComMatrix(occ; sitecolumns = sitecolumns), coords; kwargs...)

Assemblage(occ::AbstractMatrix, coords::Union{AbstractMatrix, DataFrames.DataFrame}, sites::Vector{String}, species::Vector{String}; sitecolumns = true, kwargs...) = Assemblage(ComMatrix(occ, species, sites; sitecolumns = sitecolumns), coords; kwargs...)

# a constructor that takes coords as a data.frame
function Assemblage(occ::ComMatrix, coords::DataFrames.DataFrame; kwargs...)

    if DataFrames.ncol(coords) == 2 && all(map(x -> x<:Number, eltypest(coords)))
        coords = convert(Matrix, coords)
    elseif DataFrames.ncol(coords) == 3
        xind, yind = guess_xycols(coords)
        siteind = setdiff(1:3, [xind, yind])[1]
        coords = convert(Matrix, coords[!,[xind, yind]])
    else
        error("coords must be a DataFrame with a column for sites and two columns for coordinates")
    end

    Assemblage(occ, coords; kwargs...)
end

function Assemblage(occ::ComMatrix, coords::AbstractMatrix;
      dropemptyspecies::Bool = false, dropemptysites::Bool = false, match_to_coords = true,
      traits = DataFrames.DataFrame(name = speciesnames(occ)), sitestat = DataFrames.DataFrame(sites = sitenames(occ)),
      cdtype::coordstype = auto)

    if match_to_coords
        occ, coords, sitestat = match_commat_coords(occ, coords, sitestat)
    end

    Assemblage(createLocations(coords, cdtype, sitestat), SpeciesData(occ, traits))
  end

function Assemblage(site::P, occ::SpeciesData{D};
    dropemptyspecies::Bool = false, dropemptysites::Bool = false) where {D <: Real, P <: SELocations}

    if dropemptyspecies
        dropspecies!(occ)
    end
    if dropemptysites
        dropsites!(occ, site)
    end
    Assemblage{D, P}(site, occ)
end

function createLocations(coords::AbstractMatrix, cdtype::coordstype = auto,  #by design, this is not type stable, but maybe that is OK for type constructors
        sitestat = DataFrames.DataFrame(sites = 1:size(coords, 1)))

    cdtype == pointdata && return Locations{PointData}(PointData(coords), sitestat)
    cdtype == griddata && return Locations{GridData}(GridData(coords), sitestat)
    if cdtype == auto
        try
            return Locations{GridData}(GridData(coords), sitestat)
        catch
            return Locations{PointData}(PointData(coords), sitestat)
        end
    end
end

function ComMatrix(occ::DataFrames.DataFrame; sitecolumns = true)
    if DataFrames.ncol(occ) == 3 && eltypest(occ)[3] <: AbstractString
        println("Data format identified as Phylocom")
        sites = unique(occ[!,1])
        species = unique(occ[!,3])
        js = indexin(occ[!,1], sites)
        is = indexin(occ[!,3], species)
        occ = maximum(occ[!,2]) == 1 ? sparse(is, js, true) : sparse(is, js, occ[2])
        return ComMatrix(occ, string.(collect(species)), string.(collect(sites)))
    end

    if eltypet(occ[!,1]) <: AbstractString
        species = string.(occ[:,1])
        occ = occ[!,2:end]
        sites = string.(collect(names(occ)))
    else
        species = string.(1:DataFrames.nrow(occ))
        sites = string.(collect(names(occ)))
    end

    try
        occ = dataFrametoSparseMatrix(occ, Bool)
        println("Matrix data assumed to be presence-absence")
    catch
        try
            occ = dataFrametoSparseMatrix(occ, Int) #TODO These lines mean that this code is not completely type stable.
            println("Matrix data assumed to be abundances, minimum $(minimum(occ)), maximum $(maximum(occ))")
        catch
            occ = dataFrametoSparseMatrix(occ, Float64)
            println("Matrix data assumed to be relative abundances, minimum $(minimum(occ)), maximum $(maximum(occ))")
            (minimum(occ) < 0 || maximum(occ) > 1) && @info("Values don't fall in the 0,1 range")
        end
    end

    if !sitecolumns
        return ComMatrix(collect(occ'), sites, species)
    end

    ComMatrix(occ, species, sites)
end

function ComMatrix(; species::AbstractVector = error("keyword `species` must be specified"), sites::AbstractVector = error("keyword `sites` must be specified"), abundances = 0)
    usites = unique(sites)
    uspecies = unique(species)
    js = indexin(sites, usites)
    is = indexin(species, uspecies)
    occ = abundances isa AbstractVector && length(unique(abundances)) > 1 ? sparse(is, js, abundances) : sparse(is, js, true)
    return ComMatrix(occ, string.(collect(uspecies)), string.(collect(usites)))
end

ComMatrix(occurrences::AbstractMatrix, speciesnames, sitenames; sitecolumns = true) =
    ComMatrix(occurrences; speciesnames = speciesnames, sitenames = sitenames, sitecolumns = sitecolumns)

function ComMatrix(occs; speciesnames = :auto, sitenames = :auto, sitecolumns = true)
    occurrences = sitecolumns ? occs : occs'

    if sitenames == :auto
        sitenames = ["site$i" for i in 1:size(occurrences, 2)]
    end
    if speciesnames == :auto
        speciesnames = ["species$i" for i in 1:size(occurrences, 1)]
    end

    length(speciesnames) == size(occurrences, 1) || throw(ArgumentError("length of speciesnames ($(length(speciesnames))) different from number of species in occurrences ($(size(occurrences, 1)))"))
    length(sitenames) == size(occurrences, 2) || throw(ArgumentError("length of sitenames ($(length(sitenames))) different from number of sites in occurrences ($(size(occurrences, 2)))"))
    ComMatrix{eltype(occurrences)}(sparse(occurrences), string.(speciesnames), string.(sitenames))
end

function GridData(coords::AbstractMatrix{<:Union{AbstractFloat, Missing}})
    grid = creategrid(coords)
    indices = getindices(coords, grid)
    GridData(indices, grid)
end

GridTopology(min_x, max_x, cellsize_x, min_y, max_y, cellsize_y) =
    GridTopology(range(min_x, stop = max_x, step = cellsize_x), range(min_y, stop = max_y, step = cellsize_y))


text:
module ImageNoise

using Reexport

include("NoiseAPI/NoiseAPI.jl")
include("ApplyNoise/ApplyNoise.jl")
include("ReduceNoise/ReduceNoise.jl")

@reexport using .ApplyNoise
@reexport using .ReduceNoise

end # module


text:
# ---
# title: JuliaRL\_Minimax\_OpenSpiel(tic_tac_toe)
# cover: assets/logo.svg
# description: Minimax applied to OpenSpiel(tic_tac_toe)
# date: 2021-05-22
# author: "[Jun Tian](https://github.com/findmyway)"
# ---

#+ tangle=true
using ReinforcementLearning
using OpenSpiel

function RL.Experiment(::Val{:JuliaRL}, ::Val{:Minimax}, ::Val{:OpenSpiel}, game;)
    env = OpenSpielEnv(string(game))
    agents = MultiAgentManager(
        NamedPolicy(0 => MinimaxPolicy()),
        NamedPolicy(1 => MinimaxPolicy()),
    )
    hooks = MultiAgentHook(0 => TotalRewardPerEpisode(), 1 => TotalRewardPerEpisode())
    description = "# Play `$game` in OpenSpiel with Minimax"
    Experiment(agents, env, StopAfterEpisode(1, is_show_progress=!haskey(ENV, "CI")), hooks, description)
end

#+ tangle=false
using Plots
ex = E`JuliaRL_Minimax_OpenSpiel(tic_tac_toe)`
run(ex)

# The final reward is:

ex.hook


text:
struct AxisArrayChecks{T}
    AxisArrayChecks{T}() where {T} = new{T}()
    AxisArrayChecks() = AxisArrayChecks{Union{}}()
end

struct CheckedAxisLengths end
checked_axis_lengths(::AxisArrayChecks{T}) where {T} = AxisArrayChecks{Union{T,CheckedAxisLengths}}()
check_axis_length(ks, inds, ::AxisArrayChecks{T}) where {T >: CheckedAxisLengths} = nothing
function check_axis_length(ks, inds, ::AxisArrayChecks{T}) where {T}
    if length(ks) != length(inds)
        throw(DimensionMismatch(
            "keys and indices must have same length, got length(keys) = $(length(ks))" *
            " and length(indices) = $(length(inds)).")
        )
    end
    return nothing
end

struct CheckedUniqueKeys end
checked_unique_keys(::AxisArrayChecks{T}) where {T} = AxisArrayChecks{Union{T,CheckedUniqueKeys}}()
check_unique_keys(ks, ::AxisArrayChecks{T}) where {T >: CheckedUniqueKeys} = nothing
function check_unique_keys(ks, ::AxisArrayChecks{T}) where {T}
    if allunique(ks)
        return nothing
    else
        error("All keys must be unique")
    end
end
struct CheckedOffsets end
checked_offsets(::AxisArrayChecks{T}) where {T} = AxisArrayChecks{Union{T,CheckedOffsets}}()
check_offsets(ks, inds, ::AxisArrayChecks{T}) where {T >: CheckedOffsets} = nothing
function check_offsets(ks, inds, ::AxisArrayChecks{T}) where {T}
    if firstindex(inds) === firstindex(ks)
        return nothing
    else
        throw(ArgumentError("firstindex of $ks and $inds are not the same."))
    end
end

const NoChecks = AxisArrayChecks{Union{CheckedAxisLengths,CheckedUniqueKeys,CheckedOffsets}}()


text:
# This piece of code creates a helper macro to create exceptions
# that will be thrown when the main macro is not used correctly

# Did you watch Inception, by the way?

"A dictionary (`Symbol => Any`) containing all exceptions used in this package."
EXCEPTIONS = Dict{Symbol, Any}()

"""
    @aux(
        exception_name::Symbol,
        docstring::Union{Expr, String},
        error_message_bits::Tuple{Vararg{Union{Expr, String}}},
    ) -> Expr

Create an exception with no fields and push it to the exceptions dictionary.

# Arguments
- `exception_name::Symbol`: name of the exception
- `docstring::Union{Expr, String}`: documentation string
- `error_message_bits::Tuple{Vararg{Union{Expr, String}}}`: strings and expressions which
  will be interpolated in the `showerror` output

# Returns
- `Expr`: an exception definition (struct + `showerror` overload)

See also: [`@exception`](@ref)
"""
macro aux(
    exception_name::Symbol,
    docstring::Union{Expr, String},
    error_message_bits::Union{Expr, String}...,
)
    module_name = __module__
    error_header = "$(module_name).$(exception_name):"

    return esc(
        quote
            @doc $(docstring)
            mutable struct $(exception_name) <: Exception
            end

            Base.showerror(io::IO, e::$(exception_name)) =
            print(
                io, string(
                    '\n', '\n',
                    $(error_header), '\n',
                    $(error_message_bits...), '\n',
                )
            )

            EXCEPTIONS[$(QuoteNode(exception_name))] = $(exception_name)
        end
    )
end

@aux(
    DocstringIsNotAString,
    "Exception thrown when the passed expression for docstring does not yield a string.",
    "The passed expression for docstring does not yield a string."
)

@aux(
    ErrorMessageIsNotAString,
    "Exception thrown when the passed expression(s) for error message (do)es " *
    "not yield a string.",
    "The passed expression(s) for error message (do)es not yield a string."
)

@aux(
    FieldsOnly,
    "Exception thrown when the passed argument is not a structure field.",
    "The passed argument is not a structure field (e.g. `file::String` or `file`)."
)

@aux(
    OnlyOneEquation,
    "Exception thrown when a non-context equation has been passed.",
    "Equations are only available for expressions with the first argument " *
    "equal to `context`.",
)


text:
# Autogenerated wrapper script for DASKR_jll for aarch64-linux-musl-libgfortran3
export libdaskr

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("DASKR")
JLLWrappers.@declare_library_product(libdaskr, "libdaskr.so")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libdaskr,
        "lib/libdaskr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#
# @description Unit tests for null (identity) block-cipher.
#
# @author Alejandro Claro (alejandro.claro@gmail.com)
#
# Copyright 2015 All rights reserved.
# Use of this source code is governed by a MIT-style license that can be found in the LICENSE file.
#
@testset "Null cipher tests" begin
  cipher = NullBlockCipher()

  @test block_size(cipher) == 16
  @test key_size(cipher)   == 16

  cipher = NullBlockCipher(4, 24)

  @test block_size(cipher) == 4
  @test key_size(cipher)   == 24

  @test encrypt(cipher, [ 0x01, 0x02, 0x03, 0x04 ]) == [ 0x01, 0x02, 0x03, 0x04 ]
  @test encrypt(cipher, [ 0x05, 0x06, 0x07, 0x08 ]) == [ 0x05, 0x06, 0x07, 0x08 ]

  @test_throws ArgumentError encrypt(cipher, UInt8[])
  @test_throws ArgumentError encrypt(cipher, [ 0x05, 0x06 ])
  @test_throws ArgumentError encrypt(cipher, [ 0x05, 0x06, 0x07, 0x08, 0x09 ])

  @test decrypt(cipher, [ 0x01, 0x02, 0x03, 0x04 ]) == [ 0x01, 0x02, 0x03, 0x04 ]
  @test decrypt(cipher, [ 0x05, 0x06, 0x07, 0x08 ]) == [ 0x05, 0x06, 0x07, 0x08 ]

  @test_throws ArgumentError decrypt(cipher, UInt8[])
  @test_throws ArgumentError decrypt(cipher, [ 0x05, 0x06 ])
  @test_throws ArgumentError decrypt(cipher, [ 0x05, 0x06, 0x07, 0x08, 0x09 ])
end


text:
export WingStrip, RigidLink, FuselageCut, LumpedMass, Beam, Aircraft, Airfoil, PointDrag,
	add_beam!, add_link!, add_mass!, add_wing_strip!, add_point_drag!, add_fuselage_cut!, flat_plate,
	Engine, add_engine!, add_point!, interpolate_wing!, interpolate_fuselage!,
	plain_flap, slotted_flap, slat, fixed_slat, Kruger_slat,
	get_massic_properties,
	WingSectInfo, FuselageCutInfo, interpolate!

"""
```
	struct LumpedMass{Fg <: Real}
		m::Fg
		I::Matrix{Fg}
	end
```

Lumped mass struct
"""
struct LumpedMass{Fg <: Real}
	m::Fg
	I::Matrix{Fg}
end

import Base: *
"""
Matrix multiplications for lumped mass
"""
function *(LM::LumpedMass{Fg}, M::Matrix{Fc}) where {Fg, Fc <: Real}

	LM.I .= LM.I * M

	LM

end

"""
Matrix multiplications for lumped mass
"""
function *(M::Matrix{Fc}, LM::LumpedMass{Fg}) where {Fg, Fc <: Real}

	LM.I .= M * LM.I

	LM

end

"""
Multiplication between lumped mass and scalar
"""
*(Sc::Fc, M::LumpedMass{Fg}) where {Fg, Fc <: Real} = LumpedMass{Fg}(
	m * Sc,
	M.I .* Sc
)

"""
Multiplication between lumped mass and scalar
"""
*(M::LumpedMass{Fg}, Sc::Fc) where {Fg, Fc <: Real} = Sc * M

import Base: /
"""
Division between lumped mass and scalar
"""
/(M::LumpedMass{Fg}, Sc::Fc) where {Fg, Fc <: Real} = LumpedMass{Fg}(
	M.m / Sc,
	M.I ./ Sc
)

import Base: +
"""
Addition between two lumped masses
"""
+(M1::LumpedMass{Fg}, M2::LumpedMass{Fg}) where Fg = LumpedMass{Fg}(
	M1.m + M2.m,
	M1.I .+ M2.I
)

"""
Abstract type for a fuselage cut or wing strip
"""
abstract type AbstractAeroModel end

"""
```
	struct Airfoil{Fg <: Real}
		CL0::Fg
		CLα::Fg
		Cm0::Fg
		Cmα::Fg
		polar::Vector{Fg}
		CLmin::Fg
		CLmax::Fg
	end
```

Data type defining an airfoil

* `polar` should be an array `[a, b, c]` such that `CD = CL ^ 2 * a + CL * b + c`
"""
struct Airfoil{Fg <: Real}
	CL0::Fg
	CLα::Fg
	Cm0::Fg
	Cmα::Fg
	polar::Vector{Fg}
	CLmin::Fg
	CLmax::Fg
end

Base.:+(afl1::Airfoil, afl2::Airfoil) = Airfoil(
    (afl1.CL0 + afl2.CL0),
    (afl1.CLα + afl2.CLα),
    (afl1.Cm0 + afl2.Cm0),
    (afl1.Cmα + afl2.Cmα),
    (afl1.polar .+ afl2.polar),
	(afl1.CLmin + afl2.CLmin),
	(afl1.CLmax + afl2.CLmax)
)
Base.:*(η::Number, afl::Airfoil) = Airfoil(
    η * afl.CL0,
    η * afl.CLα,
    η * afl.Cm0,
    η * afl.Cmα,
    η .* afl.polar,
	η * afl.CLmin,
	η * afl.CLmax
)
Base.:*(afl::Airfoil, η::Number) = η * afl

Airfoil(
	;
	CL0::Fg = 0.0,
	CLα::Fg = 2 * π,
	Cm0::Fg = 0.0,
	Cmα::Fg = 0.0,
	polar::Vector{Fg} = zeros(3),
	CLmin::Fg = -0.7,
	CLmax::Fg = 0.7
) where Fg = Airfoil{Fg}(
	CL0,
	CLα,
	Cm0,
	Cmα,
	polar,
	CLmin,
	CLmax
)

"""
```
	struct WingStrip{Fg <: Real} <: AbstractAeroModel
		ipt1::Int64
		ipt2::Int64
		ibeam::Int64
		Mtosys::Matrix{Fg}
		Mtouni::Matrix{Fg}
		L::Fg
		c::Fg
		xCG::Fg
		e::Fg
		afl::Airfoil{Fg}
	end
```

Wing strip data type
"""
struct WingStrip{Fg <: Real} <: AbstractAeroModel
	ipt1::Int64
	ipt2::Int64
	ibeam::Int64
	Mtosys::Matrix{Fg}
	Mtouni::Matrix{Fg}
	L::Fg
	c::Fg
	xCG::Fg
	e::Fg
	afl::Airfoil{Fg}
end

"""
```
	struct FuselageCut{Fg <: Real} <: AbstractAeroModel
		ipt1::Int64
		ipt2::Int64
		ibeam::Int64
		Mtosys::Matrix{Fg}
		Mtouni::Matrix{Fg}
		A1::Fg
		A2::Fg
		L::Fg
		Cf::Fg
	end
```

Fuselage cut data type
"""
struct FuselageCut{Fg <: Real} <: AbstractAeroModel
	ipt1::Int64
	ipt2::Int64
	ibeam::Int64
	Mtosys::Matrix{Fg}
	Mtouni::Matrix{Fg}
	A1::Fg
	A2::Fg
	L::Fg
	Cf::Fg
end

"""
```
	struct Beam{Fg <: Real}
		ipt1::Int64
		ipt2::Int64
		K::Matrix{Fg}
		g::Fg
		Mtosys::Matrix{Fg}
		Mtouni::Matrix{Fg}
	end
```

Struct to designate a simple beam element
"""
struct Beam{Fg <: Real}
	ipt1::Int64
	ipt2::Int64
	K::Matrix{Fg}
	g::Fg
	Mtosys::Matrix{Fg}
	Mtouni::Matrix{Fg}
end

"""
```
	struct PointDrag{Fg <: Real}
		ipt::Int64
		CD::Fg
		A::Fg
		V::Fg
	end
```

Struct to designate a punctual drag singularity
"""
struct PointDrag{Fg <: Real}
	ipt::Int64
	CD::Fg
	A::Fg
	V::Fg
end

"""
```
	struct Engine{Fg <: Real}
		ipt::Int64
		ξ⃗::AbstractVector
		R::Fg
	end
```

Struct to designate an engine
"""
struct Engine{Fg <: Real}
	ipt::Int64
	ξ⃗::AbstractVector
	R::Fg
end

"""
Struct to designate a perfectly rigid link
between two points
"""
struct RigidLink
	i1::Int64
	i2::Int64

	RigidLink(
		i1::Int,
		i2::Int = 0
	) = new(i1, i2)
end

"""
```
	mutable struct Aircraft{Fg <: Real}
		points::Matrix{Fg}
		wing_strips::Vector{WingStrip{Fg}}
		fuselage_cuts::Vector{FuselageCut{Fg}}
		point_drags::Vector{PointDrag{Fg}}
		engines::Vector{Engine{Fg}}
		beams::Vector{Beam{Fg}}
		masses::Vector{LumpedMass{Fg}}
		links::Vector{RigidLink}
	end
```

Mutable struct to designate an aircraft
"""
mutable struct Aircraft{Fg <: Real}
	points::Matrix{Fg}
	wing_strips::Vector{WingStrip{Fg}}
	fuselage_cuts::Vector{FuselageCut{Fg}}
	point_drags::Vector{PointDrag{Fg}}
	engines::Vector{Engine{Fg}}
	beams::Vector{Beam{Fg}}
	masses::Vector{LumpedMass{Fg}}
	links::Vector{RigidLink}
end

"""
```
	Aircraft{Fg}(points::Matrix{Fg})
```

Constructor for an aircraft, initially devoid of wings, fuselage, beams or mass points.
Takes a matrix of points (shape `(3, npts)`)
"""
Aircraft{Fg}(points::Matrix{Fg}) where Fg = Aircraft{Fg}(
	points,
	WingStrip{Fg}[],
	FuselageCut{Fg}[],
	PointDrag{Fg}[],
	Engine{Fg}[],
	Beam{Fg}[],
	[
		LumpedMass{Fg}(
			0.0,
			zeros(Fg, 3, 3)
		) for i = 1:size(points, 2)
	],
	RigidLink[]
)
Aircraft(points::Matrix{Fg}) where {Fg <: Real} = Aircraft{Fg}(points)
"""
```
	Aircraft()
```

Constructor for an aircraft initially devoid of all elements
"""
Aircraft() = Aircraft(Matrix{Float64}(undef, 3, 0))

"""
```
	function add_mass!(
		acft::Aircraft{Fg},
		ipt::Int64,
		m::Fg;
		Ixx::Fg = 0.0,
		Iyy::Fg = 0.0,
		Izz::Fg = 0.0,
		Ixy::Fg = 0.0,
		Ixz::Fg = 0.0,
		Iyz::Fg = 0.0
	) where {Fg}
```

Add point mass to an aircraft.

* `acft`: aircraft at hand
* `ipt`: index for point in which the mass will be added
* `m`: mass
* `Ixx, Iyy, Izz, Ixy, Ixz, Iyz`: moments of inertia.
	Default to zero if absent
"""
function add_mass!(
	acft::Aircraft{Fg},
	ipt::Int64,
	m::Fg;
	Ixx::Fg = 0.0,
	Iyy::Fg = 0.0,
	Izz::Fg = 0.0,
	Ixy::Fg = 0.0,
	Ixz::Fg = 0.0,
	Iyz::Fg = 0.0
) where {Fg}

	acft.masses[ipt] += LumpedMass{Fg}(
		m,
		[
			Ixx Ixy Ixz;
			Ixy Iyy Iyz;
			Ixz Iyz Izz
		]
	)

end

"""
```
	add_link!(
		acft::Aircraft{Fg},
		i1::Int,
		i2::Int = 0
	) where Fg = push!(acft.links, RigidLink(i1, i2))
```

Add rigid link to an aircraft struct
"""
add_link!(
	acft::Aircraft{Fg},
	i1::Int,
	i2::Int = 0
) where Fg = push!(acft.links, RigidLink(i1, i2))

"""
```
	function add_beam!(
		acft::Aircraft{Fg},
		ipt1::Int64,
		ipt2::Int64;
		m::Fg = 1e-3,
		ŷ::Vector{Fg} = [1.0, 0.0, 0.0],
		EA::Fg = Fmax,
		EIy::Fg = Fmax,
		EIz::Fg = Fmax,
		GJ::Fg = Fmax,
		Ixx::Fg = 1e-3,
		g::Fg = 1e-2
	) where Fg
```

Add a beam to the aircraft model

* `acft`: aircraft struct
* `ipt1`: index for first point
* `ipt2`: index for second point
* `m`: linear mass density
* `EA`: normal stiffness
* `EIy`: y axis flexural stiffness
* `EIz`: z axis flexural stiffness
* `GJ`: x axis torsional stiffness
* `Ixx`: moment of inertia around beam axis (per unit length)
* `ŷ`: y direction used to define the sectional coordinate system
* `g`: structural damping coefficient (damping matrix is set to `- Kjg`)
"""
function add_beam!(
	acft::Aircraft{Fg},
	ipt1::Int64,
	ipt2::Int64;
	m::Fg = 1e-3,
	ŷ::Vector{Fg} = [1.0, 0.0, 0.0],
	EA::Fg = Fmax,
	EIy::Fg = Fmax,
	EIz::Fg = Fmax,
	GJ::Fg = Fmax,
	Ixx::Fg = 1e-3,
	g::Fg = 1e-2
) where Fg

	x̂ = acft.points[:, ipt2] .- acft.points[:, ipt1]

	L = norm(x̂)
	x̂ ./= L

	ŷ .-= x̂ * (x̂ ⋅ ŷ)
	ny = norm(ŷ)
	if ny < 1e-5
		ŷ = [0.0, 1.0, 0.0]

		ŷ .-= x̂ * (x̂ ⋅ ŷ)
		ny = norm(ŷ)
	end
	ŷ ./= ny

	ẑ = cross(x̂, ŷ)

	Mtouni = [
		x̂ ŷ ẑ
	]
	Mtosys = inv(Mtouni)

	K = beam_get_K(
		L,
		EA,
		GJ,
		EIy,
		EIz
	)

	push!(acft.beams, Beam{Fg}(ipt1, ipt2, K, g, Mtosys, Mtouni))

	Ilat = m * L ^ 3 / 3

	LM = LumpedMass{Fg}(
		m * L,
		Mtouni * [
			(Ixx * L) 0.0 0.0;
			0.0 Ilat 0.0;
			0.0 0.0 Ilat
		] * Mtosys
	)

	acft.masses[ipt1] += LM / 2
	acft.masses[ipt2] += LM / 2

	acft.beams[end]

end

"""
```
	function add_point_drag!(
		acft::Aircraft{Fg},
		ipt::Int64,
		V::Fg,
		A::Union{Fg, Nothing} = nothing;
		CD::Fg = 2.0
	) where Fg
```

Function to add a punctual drag source to an aerodynamic model

* `acft`: aircraft model
* `ipt`: point number in aircraft mesh
* `V`: volume
* `A`: surface area (optional, assumed spherical if absent)
* `CD`: drag coeficient in relationship to `A`. Defaults to 2
	(detached flow all throughout the surface)
"""
function add_point_drag!(
	acft::Aircraft{Fg},
	ipt::Int64,
	V::Fg,
	A::Union{Fg, Nothing} = nothing;
	CD::Fg = 2.0
) where Fg

	push!(
		acft.point_drags,
		PointDrag{Fg}(
			ipt,
			CD,
			(isnothing(A) ? 4 * π * (3 * V / (π * 4)) ^ (2 / 3) : A),
			V
		)
	)

	acft.point_drags[end]

end

"""
```
	function add_engine!(
		acft::Aircraft{Fg},
		ipt::Int64,
		R::Fg;
		ξ⃗::AbstractVector = [- 1.0, 0.0, 0.0]
	) where Fg
```

Add engine to aircraft model.

* `acft`: aircraft model
* `ipt`: index for point in aircraft model
* `R`: actuator disk radius
* `ξ⃗`: vector pointing in the direction of the traction and moment vectors.
	Defaults to the opposite of the x vector
"""
function add_engine!(
	acft::Aircraft{Fg},
	ipt::Int64,
	R::Fg;
	ξ⃗::AbstractVector = [- 1.0, 0.0, 0.0]
) where Fg

	push!(
		acft.engines,
		Engine{Fg}(
			ipt,
			ξ⃗ ./ norm(ξ⃗),
			R
		)
	)

	length(acft.engines)

end

"""
Default, flat plate airfoil
"""
flat_plate = Airfoil{Float64}(0.0, 2 * π, 0.0, 0.0, [0.0, 0.0, 0.0], -0.7, 0.7)

"""
```
	function add_wing_strip!(
		acft::Aircraft{Fg},
		ipt1::Int64,
		ipt2::Int64,
		c::Fg;
		incidence::Fg = 0.0,
		afl::Airfoil{Fg} = flat_plate,
		xCG::Fg = 0.25,
		e::Fg = 0.0,
		m::Fg = 1e-3,
		ŷ::Vector{Fg} = [1.0, 0.0, 0.0],
		EA::Fg = Fmax,
		EIy::Fg = Fmax,
		EIz::Fg = Fmax,
		GJ::Fg = Fmax,
		Ixx::Fg = 1e-3,
		g::Fg = 1e-2
	) where Fg
```

Add wing strip to aircraft model

* `acft`: aircraft struct
* `ipt1`: index for point number 1
* `ipt2`: index for point number 2
* `c`: chord
* `incidence`: incidence (degrees)
* `afl`: airfoil
* `xCG`: position of the center of gravity (perc. of chord from aerodynamic center)
* `e`: position of the elastic center (perc. of chord from aerodynamic center)
* `m`: linear density of the wing section
* `ŷ`: y vector, perpendicular to the axis of the wing,
with which the flexural resistances `EIy` and `EIz` are calculated.
Defaults to x axis
* `EA`: axial resistance for the beam
* `EIy`: flexural resistance in the axis indicated by `ŷ`
* `EIz`: flexural resistance in the axis perpendicular to the wing and
the `ŷ` axis
* `GJ`: torsional resistance
* `Ixx`: inertia of the wing section around its elastic center (per unit length)
* `g`: damping coefficient for section beam
"""
function add_wing_strip!(
	acft::Aircraft{Fg},
	ipt1::Int64,
	ipt2::Int64,
	c::Fg;
	incidence::Fg = 0.0,
	afl::Airfoil{Fg} = flat_plate,
	xCG::Fg = 0.25,
	e::Fg = 0.0,
	m::Fg = 1e-3,
	ŷ::Vector{Fg} = [1.0, 0.0, 0.0],
	EA::Fg = Fmax,
	EIy::Fg = Fmax,
	EIz::Fg = Fmax,
	GJ::Fg = Fmax,
	Ixx::Fg = 1e-3,
	g::Fg = 1e-2
) where Fg

	_afl = Airfoil{Fg}(
		afl.CL0 + afl.CLα * deg2rad(incidence),
		afl.CLα,
		afl.Cm0 + afl.Cmα * deg2rad(incidence),
		afl.Cmα,
		copy(afl.polar),
		afl.CLmin,
		afl.CLmax
	)

	b = add_beam!(
		acft,
		ipt1,
		ipt2;
		m = m,
		Ixx = Ixx,
		ŷ = ŷ,
		EA = EA,
		EIy = EIy,
		EIz = EIz,
		GJ = GJ,
		g = g
	)

	sect = WingStrip{Fg}(
		ipt1, ipt2, length(acft.beams),
		b.Mtosys,
		b.Mtouni,
		norm(acft.points[:, ipt2] .- acft.points[:, ipt1]),
		c,
		xCG,
		e,
		_afl
	)

	push!(acft.wing_strips, sect)

	sect

end

"""
```
	function add_fuselage_cut!(
		acft::Aircraft{Fg},
		ipt1::Int64,
		ipt2::Int64,
		A1::Fg,
		A2::Fg;
		m::Fg = 1e-3,
		ŷ::Vector{Fg} = [0.0, 1.0, 0.0],
		EA::Fg = Fmax,
		EIy::Fg = Fmax,
		EIz::Fg = Fmax,
		GJ::Fg = Fmax,
		Ixx::Fg = 1e-3,
		Cf::Fg = 0.0,
		g::Fg = 1e-2
	) where Fg
```

Add fuselage section to aerodynamic model

* `acft`: aircraft struct
* `ipt1`: index for point number 1
* `ipt2`: index for point number 2
* `A1`: sectional area at point 1
* `A2`: sectional area at point 2
* `m`: linear density of the wing section
* `ŷ`: y vector, perpendicular to the axis of the fuselage,
with which the flexural resistances `EIy` and `EIz` are calculated.
Defaults to y axis
* `EA`: axial resistance for the beam
* `EIy`: flexural resistance in the axis indicated by `ŷ`
* `EIz`: flexural resistance in the axis perpendicular to the fuselage and
the `ŷ` axis
* `GJ`: torsional resistance
* `Ixx`: inertia of the fuselage section around its axis (per unit length)
* `Cf`: friction coefficient
* `g`: damping coefficient for associated beam element
"""
function add_fuselage_cut!(
	acft::Aircraft{Fg},
	ipt1::Int64,
	ipt2::Int64,
	A1::Fg,
	A2::Fg;
	m::Fg = 1e-3,
	ŷ::Vector{Fg} = [0.0, 1.0, 0.0],
	EA::Fg = Fmax,
	EIy::Fg = Fmax,
	EIz::Fg = Fmax,
	GJ::Fg = Fmax,
	Ixx::Fg = 1e-3,
	Cf::Fg = 0.0,
	g::Fg = 1e-2
) where Fg

	b = add_beam!(
		acft,
		ipt1,
		ipt2;
		m = m,
		ŷ = ŷ,
		EA = EA,
		EIy = EIy,
		EIz = EIz,
		GJ = GJ,
		Ixx = Ixx,
		g = g
	)

	L = norm(acft.points[:, ipt1] .- acft.points[:, ipt2])

	cut = FuselageCut{Fg}(
		ipt1,
		ipt2,
		length(acft.beams),
		b.Mtosys,
		b.Mtouni,
		A1,
		A2,
		L,
		Cf
	)

	push!(acft.fuselage_cuts, cut)

	cut

end

"""
```
	function add_point!(
		acft::Aircraft{Fg},
		point::AbstractVector,
		tol::Real = 1e-10
	) where Fg
```

Add a point to an aircraft.

* `acft`: aircraft at hand
* `point`: point at hand (vector of length 3)
* `tol`: tolerance for point merging
"""
function add_point!(
	acft::Aircraft{Fg},
	point::AbstractVector,
	tol::Real = 1e-10
) where Fg

	if tol != 0.0 && size(acft.points, 2) != 0
		indmin, distmin = let dp = acft.points .- point
			dists = sqrt.(dp[1, :] .^ 2 .+ dp[2, :] .^ 2 .+ dp[3, :] .^ 2)

			indmin = argmin(dists)

			indmin, dists[indmin]
		end

		if distmin < tol
			return indmin
		end
	end

	acft.points = [acft.points point]
	push!(acft.masses, LumpedMass{Fg}(zero(Fg), zeros(Fg, 3, 3)))

	size(acft.points, 2)

end

"""
```
	function interpolate_wing!(
		acft::Aircraft,
		dist::AbstractVector,
		point1::AbstractVector,
		point2::AbstractVector,
		tol::Real = 1e-10;
		args1 = [],
		args2 = [],
		kwargs1 = Dict(),
		kwargs2 = Dict()
	)
```

Add a wing by interpolating properties between two points

* `acft`: the aircraft at hand
* `dist`: array with points in a 0 to 1 interval, used
	as base for interpolation
* `point1`: coordinates for first point
* `point2`: coordinates for second point
* `tol`: optional argument for point merging tolerance
* `args1`: arguments for first point
* `args2`: arguments for second point
* `kwargs1`: keyword arguments for first point
* `kwargs2`: keyword arguments for second point

Returns: array of point indices and array of wing section indices
"""
function interpolate_wing!(
	acft::Aircraft,
	dist::AbstractVector,
	point1::AbstractVector,
	point2::AbstractVector,
	tol::Real = 1e-10;
	args1 = [],
	args2 = [],
	kwargs1 = Dict(),
	kwargs2 = Dict()
)

	point_inds = [
		add_point!(
			acft,
			(@. point1 * (1.0 - η) + point2 * η),
			tol
		) for η in dist
	]

	distmed = (dist[1:(end - 1)] .+ dist[2:end]) ./ 2

	sect_inds = [
		begin
			args = [
				(a1 * (1.0 - η) + a2 * η) for (a1, a2) in zip(args1, args2)
			]
			kwargs = [
				k1 => (v1 * (1.0 - η) + kwargs2[k1] * η) for (
					k1, v1
				) in kwargs1
			]

			add_wing_strip!(
				acft,
				point_inds[i],
				point_inds[i + 1],
				args...;
				kwargs...
			)

			length(acft.wing_strips)
		end for (i, η) in enumerate(distmed)
	]

	point_inds, sect_inds

end

Base.length(afl::Airfoil) = 1
Base.iterate(afl::Airfoil) = (afl, false)
Base.iterate(afl::Airfoil, state::Bool) = nothing

"""
```
	function interpolate_fuselage!(
		acft::Aircraft,
		dist::AbstractVector,
		point1::AbstractVector,
		point2::AbstractVector,
		A1::Real,
		A2::Real,
		tol::Real = 1e-10;
		args1 = [],
		args2 = [],
		kwargs1 = Dict(),
		kwargs2 = Dict()
	)
```

Add a fuselage by interpolating properties between two points

* `acft`: the aircraft at hand
* `dist`: array with points in a 0 to 1 interval, used
	as base for interpolation
* `point1`: coordinates for first point
* `point2`: coordinates for second point
* `A1`: cross-sectional area at first section
* `A2`: cross-sectional area at second section
* `tol`: optional argument for point merging tolerance
* `args1`: arguments for first point
* `args2`: arguments for second point
* `kwargs1`: keyword arguments for first point
* `kwargs2`: keyword arguments for second point

Returns: array of point indices and array of wing section indices
"""
function interpolate_fuselage!(
	acft::Aircraft,
	dist::AbstractVector,
	point1::AbstractVector,
	point2::AbstractVector,
	A1::Real,
	A2::Real,
	tol::Real = 1e-10;
	args1 = [],
	args2 = [],
	kwargs1 = Dict(),
	kwargs2 = Dict()
)

	point_inds = [
		add_point!(
			acft,
			(@. point1 * (1.0 - η) + point2 * η),
			tol
		) for η in dist
	]

	distmed = (dist[1:(end - 1)] .+ dist[2:end]) ./ 2

	sect_inds = [
		begin
			args = [
				(@. a1 * (1.0 - η) + a2 * η) for (a1, a2) in zip(args1, args2)
			]
			kwargs = [
				k1 => (@. v1 * (1.0 - η) + kwargs2[k1] * η) for (
					k1, v1
				) in kwargs1
			]

			add_fuselage_cut!(
				acft,
				point_inds[i],
				point_inds[i + 1],
				A1 * (1.0 - dist[i]) + A2 * dist[i],
				A1 * (1.0 - dist[i + 1]) + A2 * dist[i + 1],
				args...;
				kwargs...
			)

			length(acft.fuselage_cuts)
		end for (i, η) in enumerate(distmed)
	]

	point_inds, sect_inds

end

const _default_flap_deflection = 0.698
const _default_slat_deflection = 0.698

"""
```
	function plain_flap(
		afl::Airfoil{Fg},
		cperc::Real,
		deflection::Real;
		sweep::Real = 0.0
	) where Fg
```

Convert an airfoil to its approximate, plain flap equivalent,
according to Raymer's "Aircraft Design: a Conceptual Approach"

* `afl`: original airfoil
* `cperc`: percentage of chord
* `deflection`: deflection in degrees
* `sweep`: sweep in degrees (optional)
"""
function plain_flap(
	afl::Airfoil{Fg},
	cperc::Real,
	deflection::Real;
	sweep::Real = 0.0
) where Fg

	deflection = deg2rad(deflection)

	Δα = cperc * deflection

	dCLm = 0.9 * deflection * cosd(sweep) / _default_flap_deflection

	Airfoil{Fg}(
		afl.CL0 + Δα * afl.CLα,
		afl.CLα,
		afl.Cm0 + Δα * afl.Cmα,
		afl.Cmα,
		[
			afl.polar[1],
			afl.polar[2],
			afl.polar[3] + rad2deg(deflection) * 0.0023
		],
		afl.CLmin + dCLm,
		afl.CLmax + dCLm
	)

end

"""
```
	function slotted_flap(
		afl::Airfoil{Fg},
		cperc::Real,
		deflection::Real;
		sweep::Real = 0.0,
		n_slots::Int64 = 1
	) where Fg
```

Convert an airfoil to its approximate, slotted flap equivalent,
according to Raymer's "Aircraft Design: a Conceptual Approach"

* `afl`: original airfoil
* `cperc`: percentage of chord
* `deflection`: deflection in degrees
* `sweep`: sweep in degrees (optional)
* `n_slots`: number of slots. Ideally between 1 and 3. Defaults
	to 1
"""
function slotted_flap(
	afl::Airfoil{Fg},
	cperc::Real,
	deflection::Real;
	sweep::Real = 0.0,
	n_slots::Int64 = 1
) where Fg

	deflection = deg2rad(deflection)

	Δα = cperc * deflection

	mult = 1.0 + 0.3 * n_slots

	dCLm = mult * deflection * cosd(sweep) * cperc / _default_flap_deflection

	Airfoil{Fg}(
		afl.CL0 + Δα * afl.CLα,
		afl.CLα,
		afl.Cm0 + Δα * afl.Cmα,
		afl.Cmα,
		[
			afl.polar[1],
			afl.polar[2],
			afl.polar[3] + rad2deg(deflection) * 0.0023
		],
		afl.CLmin + dCLm,
		afl.CLmax + dCLm
	)

end

"""
```
	function fixed_slat(
		afl::Airfoil{Fg},
		cperc::Real,
		deflection::Real;
		sweep::Real = 0.0
	) where Fg
```

Convert an airfoil to its approximate, fixed slat equivalent,
according to Raymer's "Aircraft Design: a Conceptual Approach"

* `afl`: original airfoil
* `cperc`: percentage of chord
* `deflection`: deflection in degrees
* `sweep`: sweep in degrees (optional)
"""
function fixed_slat(
	afl::Airfoil{Fg},
	cperc::Real,
	deflection::Real;
	sweep::Real = 0.0
) where Fg

	deflection = deg2rad(deflection)

	Δα = 0.0 # cperc * deflection

	dCLm = 0.2 * deflection * cosd(sweep) / _default_slat_deflection

	Airfoil{Fg}(
		afl.CL0 + Δα * afl.CLα,
		afl.CLα,
		afl.Cm0 + Δα * afl.Cmα,
		afl.Cmα,
		[
			afl.polar[1],
			afl.polar[2],
			afl.polar[3] + rad2deg(deflection) * 0.0023
		],
		afl.CLmin + dCLm,
		afl.CLmax + dCLm
	)

end

"""
```
	function slat(
		afl::Airfoil{Fg},
		cperc::Real,
		deflection::Real;
		sweep::Real = 0.0
	) where Fg
```

Convert an airfoil to its approximate, slat equivalent,
according to Raymer's "Aircraft Design: a Conceptual Approach"

* `afl`: original airfoil
* `cperc`: percentage of chord
* `deflection`: deflection in degrees
* `sweep`: sweep in degrees (optional)
"""
function slat(
	afl::Airfoil{Fg},
	cperc::Real,
	deflection::Real;
	sweep::Real = 0.0
) where Fg

	deflection = deg2rad(deflection)

	Δα = 0.0 # cperc * deflection

	dCLm = 0.4 * cperc * deflection * cosd(sweep) / _default_slat_deflection

	Airfoil{Fg}(
		afl.CL0 + Δα * afl.CLα,
		afl.CLα,
		afl.Cm0 + Δα * afl.Cmα,
		afl.Cmα,
		[
			afl.polar[1],
			afl.polar[2],
			afl.polar[3] + rad2deg(deflection) * 0.0023
		],
		afl.CLmin + dCLm,
		afl.CLmax + dCLm
	)

end

"""
```
	function Kruger_slat(
		afl::Airfoil{Fg},
		cperc::Real,
		deflection::Real;
		sweep::Real = 0.0
	) where Fg
```

Convert an airfoil to its approximate, Kruger slat equivalent,
according to Raymer's "Aircraft Design: a Conceptual Approach"

* `afl`: original airfoil
* `cperc`: percentage of chord
* `deflection`: deflection in degrees
* `sweep`: sweep in degrees (optional)
"""
function Kruger_slat(
	afl::Airfoil{Fg},
	cperc::Real,
	deflection::Real;
	sweep::Real = 0.0
) where Fg

	deflection = deg2rad(deflection)

	Δα = 0.0 # cperc * deflection

	dCLm = 0.3 * deflection * cosd(sweep) / _default_slat_deflection

	Airfoil{Fg}(
		afl.CL0 + Δα * afl.CLα,
		afl.CLα,
		afl.Cm0 + Δα * afl.Cmα,
		afl.Cmα,
		[
			afl.polar[1],
			afl.polar[2],
			afl.polar[3] + rad2deg(deflection) * 0.0023
		],
		afl.CLmin + dCLm,
		afl.CLmax + dCLm
	)

end

"""
```
	function get_massic_properties(
		acft::Aircraft{Fg}
	) where Fg
```

Get mass properties for an aircraft.

Returns in format:

```
(
	I = [
		1.0 0.0 0.0;
		0.0 1.0 0.0;
		0.0 0.0 1.0 # example
	],
	CG = [0.0, 0.0, 0.0],
	M = 1.0
)
```
"""
function get_massic_properties(
	acft::Aircraft{Fg}
) where Fg

	I = zeros(3, 3)
	CG = zeros(3)
	M = 0.0

	for (i, m) in enumerate(acft.masses)
		r = acft.points[:, i]

		M += m.m
		I += m.I + m.m * [
			(r[2] ^ 2 + r[3] ^ 2) (r[1] * r[2]) (r[1] * r[3]);
			(r[1] * r[2]) (r[1] ^ 2 + r[3] ^ 2) (r[2] * r[3]);
			(r[1] * r[3]) (r[2] * r[3]) (r[1] ^ 2 + r[2] ^ 2)
		]
		CG += m.m .* r
	end

	(
		M = M,
		CG = CG / M,
		I = I
	)

end

struct WingSectInfo
	pt::AbstractVector
	args
	kwargs::Dict{Symbol, Any}
end

"""
```
	WingSectInfo(
		pt::AbstractVector,
		c::Real,
		args...;
		kwargs...
	)
```

Struct to contain wing section information.

* `pt`: vector with quarter chord location
* `c`: local chord

Pass any other arguments and keyword arguments to be passed
to `interpolate_wing!` when `interpolate!` is called.
"""
WingSectInfo(
	pt::AbstractVector,
	args...;
	kwargs...
) = WingSectInfo(
	pt,
	args,
	Dict{Symbol, Any}(
		kwargs...
	)
)

"""
```
	interpolate!(
		acft,
		distribution,
		s1::WingSectInfo,
		s2::WingSectInfo,
		tol::Real = 1e-10
	)
```

Interpolate two wing sections, as defined by `WingSectInfo`.

Returns the same as `interpolate_wing!`
"""
interpolate!(
	acft::Aircraft,
	distribution,
	s1::WingSectInfo,
	s2::WingSectInfo,
	tol::Real = 1e-10
) = interpolate_wing!(
	acft,
	(
		isa(distribution, Int) ?
		collect(LinRange(0.0, 1.0, distribution)) :
		collect(distribution)
	),
	s1.pt,
	s2.pt,
	tol;
	args1 = s1.args,
	args2 = s2.args,
	kwargs1 = s1.kwargs,
	kwargs2 = s2.kwargs
)

struct FuselageCutInfo
	pt::AbstractVector
	A::Real
	args
	kwargs::Dict{Symbol, Any}
end

"""
```
	FuselageCutInfo(
		pt::AbstractVector,
		A::Real,
		args...;
		kwargs...
	)
```

Struct to contain fuselage cut information.

* `pt`: vector with center location
* `A`: cross-sectional area

Pass any other arguments and keyword arguments to be passed
to `interpolate_fuselage!` when `interpolate!` is called.
"""
FuselageCutInfo(
	pt::AbstractVector,
	A::Real,
	args...;
	kwargs...
) = FuselageCutInfo(
	pt,
	A,
	args,
	Dict{Symbol, Any}(
		kwargs...
	)
)

"""
```
	interpolate!(
		acft,
		distribution,
		s1::FuselageCutInfo,
		s2::FuselageCutInfo,
		tol::Real = 1e-10
	)
```

Interpolate two wing sections, as defined by `FuselageCutInfo`.

Returns the same as `interpolate_fuselage!`
"""
interpolate!(
	acft::Aircraft,
	distribution,
	s1::FuselageCutInfo,
	s2::FuselageCutInfo,
	tol::Real = 1e-10
) = interpolate_fuselage!(
	acft,
	(
		isa(distribution, Int) ?
		collect(LinRange(0.0, 1.0, distribution)) :
		collect(distribution)
	),
	s1.pt,
	s2.pt,
	s1.A,
	s2.A,
	tol;
	args1 = s1.args,
	args2 = s2.args,
	kwargs1 = s1.kwargs,
	kwargs2 = s2.kwargs
)

"""
Interpolate property between adjacent wing sections so as to obtain
values at nodes
"""
function _wstrip_interp(
	acft::Aircraft{Fg},
	prop::Symbol
) where Fg

	cnt = zeros(Fg, size(acft.points, 2))
	v = zeros(Fg, size(acft.points, 2))

	for wstrip in acft.wing_strips
		p = getproperty(wstrip, prop)
		L = wstrip.L

		w = 1.0 / (L + 1e-10)

		cnt[wstrip.ipt1] += w
		v[wstrip.ipt1] += p * w

		cnt[wstrip.ipt2] += w
		v[wstrip.ipt2] += p * w
	end

	@. v / (cnt + 1e-10)

end

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
"""
    assemble_mass_matrix(elements, volumes[, weights])

Assemble 3D mass matrix using P1 finite elements.

This function is based on the Matlab function `mass_matrixP1_3D.m` from the Matlab
nodal matrix assembly toolbox by Jan Valdman and Talal Rahmnan:

https://uk.mathworks.com/matlabcentral/fileexchange/27826-fast-fem-assembly-nodal-elements

Talal Rahman and Jan Valdman: Fast MATLAB assembly of FEM matrices in 2D and 3D: nodal
elements, Applied Mathematics and Computation 219, 7151–7158 (2013).
"""
function assemble_mass_matrix(elements, volumes, weights = 1)
    x = kron(ones(1, 4), elements)
    y = kron(elements, ones(1, 4))
    if length(weights) == 1 || length(weights) == size(elements, 1)
        # P0 weights
        z = kron(weights .* volumes, reshape((ones(4, 4) + I(4)) / 20, 1, 16))
    else
        # P1 weights
        M1 = [6 2 2 2; 2 2 1 1; 2 1 2 1; 2 1 1 2] / 120
        M2 = M1[[4, 1, 2, 3], [4, 1, 2, 3]]
        M3 = M2[[4, 1, 2, 3], [4, 1, 2, 3]]
        M4 = M3[[4, 1, 2, 3], [4, 1, 2, 3]]
        z = (
            kron(volumes .* weights[elements[:, 1]], reshape(M1, 1, 16)) +
            kron(volumes .* weights[elements[:, 2]], reshape(M2, 1, 16)) +
            kron(volumes .* weights[elements[:, 3]], reshape(M3, 1, 16)) +
            kron(volumes .* weights[elements[:, 4]], reshape(M4, 1, 16))
        )
    end

    # Enforce symmetry
    sym(x) = (x + x') / 2

    sym(sparse(x[:], y[:], z[:]))
end


text:
"""
    ProblemParSetter(state_names,par_names,popt_names)
    ProblemParSetter(sys::ODESystem, popt_names; strip=false)

Helps keeping track of a subset of initial states and paraemters to be optimized.

# Arguments
- `state_names`: ComponentVector or Axis of all the initial states of the problem
- `par_names`: all the parameters of the problem
- `popt_names`: the parameter/initial states to be optimized.

If all of `state_names`, `par_names`, and `popt_names` are type-inferred Axes,
then also the constructed ProblemParSetter is type-inferred.

The states and parameters can be extracted from an `ModelingToolkit.ODESystem`.
If `strip=true`, then namespaces of parameteres of a composed system are removed,
e.g. `subcomp₊p` becomes `p`.
"""
struct ProblemParSetter{NOPT, POPTA <: AbstractAxis, SA <: AbstractAxis, PA <: AbstractAxis} <: AbstractProblemParSetter
    # u0_opt::NTuple{NS, Symbol}
    # p_opt::NTuple{NP, Symbol}
    ax_paropt::POPTA
    ax_state::SA
    ax_par::PA
    is_state::NTuple{NOPT, Bool}
    is_p::NTuple{NOPT, Bool}
end

function ProblemParSetter(ax_state::AbstractAxis, ax_par::AbstractAxis, ax_paropt::AbstractAxis)
    popt_names = keys(CA.indexmap(ax_paropt))
    state_names = keys(CA.indexmap(ax_state))
    par_names = keys(CA.indexmap(ax_par))
    NOPT = length(CA.indexmap(ax_paropt))
    is_state = ntuple(i -> popt_names[i] ∈ state_names, NOPT)
    is_p = ntuple(i -> popt_names[i] ∈ par_names, NOPT)
    is_u0_or_p = is_state .| is_p
    all(is_u0_or_p) || @warn(
        "missing optimization parameters in system: " * join(popt_names[collect(.!is_u0_or_p)],", "))
    is_u0_and_p = is_state .& is_p
    any(is_u0_and_p) && @warn(
        "expted parameter names and state names to be distinct, but occured in both: " * join(popt_names[collect(is_u0_and_p)],", "))
    ProblemParSetter{NOPT, typeof(ax_paropt),typeof(ax_state), typeof(ax_par)}(ax_paropt, ax_state, ax_par, is_state, is_p)
end

# function ProblemParSetter(state_names::NTuple{NS, Symbol}, par_names::NTuple{NP, Symbol}, popt_names::NTuple{NOPT, Symbol}) where {NS, NP, NOPT}
#     # not type stable
#     ProblemParSetter(Axis(state_names), Axis(par_names), Axis(popt_names))
# end

function ProblemParSetter(state_template,par_template,popt_template)
    # not type-stable
    ProblemParSetter(
        _get_axis(state_template),
        _get_axis(par_template),
        _get_axis(popt_template),
    )
end

# depr: need a full-fledged axis
# function _get_axis(x::Union{Tuple, AbstractArray})
#     Axis(Tuple(i for i in symbol.(x)))
# end
_get_axis(x::ComponentVector) = first(getaxes(x))
_get_axis(x::AbstractAxis) = x

function ProblemParSetter(sys::ODESystem,popt_names; strip=false)
    ft = strip ? strip_namespace : identity
    state_names = ft.(symbol.(states(sys)))
    par_names = ft.(symbol.(parameters(sys)))
    ProblemParSetter(Axis(state_names), Axis(par_names), _get_axis(popt_names))
end

# count_state(::ProblemParSetter{N, POPTA, SA, PA}) where {N, POPTA, SA, PA} = length(CA.indexmap(SA))
# count_par(::ProblemParSetter{N, POPTA, SA, PA}) where {N, POPTA, SA, PA} = length(CA.indexmap(PA))
# count_paropt(::ProblemParSetter{N, POPTA, SA, PA}) where {N, POPTA, SA, PA} = N

# todo change to length(ax) when this becomes availale in ComponentArrays
count_state(pset::ProblemParSetter) = lastindex(pset.ax_state) - firstindex(pset.ax_state) + 1
count_par(pset::ProblemParSetter) = lastindex(pset.ax_par) - firstindex(pset.ax_par) + 1
count_paropt(pset::ProblemParSetter) = lastindex(pset.ax_paropt) - firstindex(pset.ax_paropt) + 1

# axis_state(::ProblemParSetter{N, POPTA, SA, PA}) where {N, POPTA, SA, PA} = SA
# axis_par(::ProblemParSetter{N, POPTA, SA, PA}) where {N, POPTA, SA, PA} = PA
# axis_paropt(::ProblemParSetter{N, POPTA, SA, PA}) where {N, POPTA, SA, PA} = POPTA

axis_state(ps::ProblemParSetter) = ps.ax_state
axis_par(ps::ProblemParSetter) = ps.ax_par
axis_paropt(ps::ProblemParSetter) = ps.ax_paropt

# function _ax_symbols(ax::Union{AbstractAxis, CA.CombinedAxis}; prefix="₊")
#     # strip the first prefix, convert to symbol and retun generator
#     (i for i in _ax_string_prefixed(ax; prefix) .|> (x -> x[(sizeof(prefix)+1):end]) .|> Symbol)
# end
# function _ax_symbols_tuple(ax::Union{AbstractAxis, CA.CombinedAxis}; kwargs...)
#     Tuple(_ax_symbols(ax; kwargs...))::NTuple{CA.last_index(ax), Symbol}
# end
# function _ax_symbols_vector(ax::Union{AbstractAxis, CA.CombinedAxis}; kwargs...)
#     # strip the first prefix, convert to symbol and collect into tuple
#     collect(_ax_symbols(ax; kwargs...))::Vector{Symbol}
# end

symbols_state(pset::ProblemParSetter) = _ax_symbols_tuple(axis_state(pset))
symbols_par(pset::ProblemParSetter) = _ax_symbols_tuple(axis_par(pset))
symbols_paropt(pset::ProblemParSetter) = _ax_symbols_tuple(axis_paropt(pset))

# Using unexported interface of ComponentArrays.axis, one place to change
"Accessor function for index from ComponentIndex"
idx(ci::CA.ComponentIndex) = ci.idx

"""
    prob_new = update_statepar(pset::ProblemParSetter, popt, prob::ODEProblem)
    u0new, pnew = update_statepar(pset::ProblemParSetter, popt, u0, p)

Return an updated problem or updates states and parameters where
values corresponding to positions in `popt` are set.
"""
# function update_statepar(pset::ProblemParSetter, popt, prob::ODEProblem)
#     u0,p = update_statepar(pset, popt, prob.u0, prob.p)
#     remake(prob; u0, p)
# end

function update_statepar(pset::ProblemParSetter, popt::TO, u0::TU, p::TP) where {TO, TU, TP}
    popt_state, popt_p = _separate_state_p(pset, popt)
    u0new = length(popt_state) == 0 ? u0 : begin
        u0_l = label_state(pset, u0)
        fdata = (u0 isa ComponentVector) ? identity : getdata # stip labels?
        TUP = typeof(similar(u0_l, promote_type(eltype(TU), eltype(TO))))
        u0new = fdata(_update_cv_top(u0_l, popt_state)::TUP)
    end
    pnew = length(popt_p) == 0 ? p : begin
        p_l = label_par(pset,p)
        fdata = (p isa ComponentVector) ? identity : getdata # stip labels?
        TPP = typeof(similar(p_l, promote_type(eltype(TP), eltype(TO))))
        pnew = fdata(_update_cv_top(p_l, popt_p)::TPP)
    end
    u0new, pnew
end

# extract p and state components of popt into separate ComponentVectors
function _separate_state_p(pset, popt)
    popt_l = label_paropt(pset, popt)
    syms_p = Tuple(k for (i,k) in enumerate(keys(popt_l)) if pset.is_p[i])
    syms_s = Tuple(k for (i,k) in enumerate(keys(popt_l)) if pset.is_state[i])
    # popt_state = _get_index_axis(popt_l, _get_axis_of_lengths_for_syms(popt_l,
    #     Tuple(k for (i,k) in enumerate(keys(popt_l)) if pset.is_state[i])))
    # #popt_p = popt_l[Axis( # WAIT use proper indexing when supported with ComponentArrays
    # popt_p = _get_index_axis(popt_l, _get_axis_of_lengths_for_syms(popt_l,
    #         Tuple(k for (i,k) in enumerate(keys(popt_l)) if pset.is_p[i])))
    #popt_state, popt_p
    _indexof(popt_l, syms_s), _indexof(popt_l, syms_p)
end

# TODO after cv[()] works with empty and mono-Tuple
get_empty(cv::ComponentVector) = ComponentVector(similar(getdata(cv),0), (ComponentArrays.NullAxis(),))
function _indexof(cv::ComponentVector{T}, syms::NTuple{N,Symbol}) where {T,N}
    N == 1 && return(cv[KeepIndex(syms[1])])::ComponentVector{T}
    N == 0 && return(get_empty(cv))::ComponentVector{T}
    res = cv[syms]::ComponentVector{T}
end

# function typed_from_generator(type::Type, v0, vgen) where T
#     if type <: AbstractArray
#     end
#     typeof(v0)(vgen)
# end
# # AbstractVector{T} is not more specific than ::Type, need to support all used concrete
# #typed_from_generator(::Type{AbstractVector{T}}, v0, vgen) where T = convert(typeof(v0), v0 .+ vgen)::typeof(v0)
# typed_from_generator(::Type{Vector{T}}, v0, vgen) where T = convert(typeof(v0), v0 .+ vgen)::typeof(v0)

# typed_from_generator(v0, vgen) = typeof(v0)(vgen)
# function typed_from_generator(v0::AbstractVector, vgen)
#     # convert to std vector, because typeof(v0) does not contain all entries
#     T = Vector{eltype(v0)}
#     convert(T, v0 .+ vgen)::T
# end

# type piracy - try to get into CompponentArrays
# Base.getindex(cv::ComponentVector, ax::AbstractAxis) = _get_index_axis(cv,ax)
# Base.getindex(cv::ComponentVector, cv_template::ComponentVector) = _get_index_axis(
#     cv,first(getaxes(cv_template)))

# function get_paropt_labeled(pset::ProblemParSetter, u0, p)
#     ax = axis_paropt(pset)
#     keys_ax = keys(ax)
#     u0l = label_state(pset, u0)
#     pl = label_par(pset, p)
#     (i,k) = first(enumerate(keys_ax))
#     #(i,k) = (2, keys_ax[2])
#     fik = (i,k) -> begin
#         cvs = pset.is_state[i] ? getproperty(u0l,k) : (pset.is_p[i] ? getproperty(pl,k) : missing)
#         ismissing(cvs) && return(
#             ComponentVector(NamedTuple{(k,)}((fill(missing,length(ax[k].idx)),))))
#         !(cvs isa ComponentVector) && return(
#             cvs = pset.is_state[i] ? u0l[KeepIndex(k)] : pl[KeepIndex(k)])
#         axs = ax[k].ax
#         #@show cvs, axs, typeof(cvs)
#         # TODO replace by cvs[axs] when _get_index was merged
#         _get_index_axis(cvs, axs)
#     end
#     tmp = (fik(i,k) for (i,k) in enumerate(keys_ax))
#     T = promote_type(eltype(u0), eltype(p))
#     res = reduce(vcat, tmp)::ComponentVector{T, Vector{T}}
#     # @infiltrate
#     label_paropt(pset, res) # reattach axis for type inference
# end

function get_paropt_labeled(pset::ProblemParSetter, u0, p)
    let u0_l = label_state(pset, u0), p_l = label_par(pset, p)
        fik = (i,k) -> pset.is_state[i] ? u0_l[KeepIndex(k)] : (pset.is_p[i] ? p_l[KeepIndex(k)] : missing)
        tmp = (fik(ik[1], ik[2]) for ik in enumerate(keys(axis_paropt(pset))))
        T = promote_type(eltype(u0), eltype(p))
        res = reduce(vcat, tmp)::ComponentVector{T, Vector{T}}
        label_paropt(pset, res) # reattach axis for type inference
    end
end

# attach type in
# type piracy I - until get this into ComponentArrays
@inline CA.getdata(x::ComponentArray{T,N,A}) where {T,N,A} = getfield(x, :data)::A
@inline CA.getdata(x::ComponentVector{T,A}) where {T,A} = getfield(x, :data)::A

# # extends Base.merge to work on SVector
# # ?type piracy
# merge(x::T, y::NamedTuple) where T<:SLArray = T(merge(NamedTuple(x),y)...)

# # and on Labelled Arrays
# function merge(x::T, y::NamedTuple) where T<:LArray
#     xnew = deepcopy(x)
#     for (key,val) in pairs(y)
#         xnew[key] = val
#     end
#     xnew
# end


text:
using Plots
gr(xlims=[0,2*pi],ylims=[-1.2,1.2],border=:none)

Np_F = 4
h = 2*pi/Np_F

fourier(x) = @. cos(pi*x)
r_plot = 0:4/200:2
plot(r_plot*pi,fourier(r_plot),lw=3,legend=nothing)
plot!(r_plot*pi,fourier(r_plot,2),lw=3,legend=nothing)
plot!(r_plot*pi,fourier(r_plot,4),lw=3,legend=nothing)
# plot!(r_plot*pi,fourier(r_plot,8),lw=3,legend=nothing)

# y = ones(size(r_plot))
# l = 0.3
# u = 0.3
# plot!([-pi,0],[2,2],fill=(0,:grey),fillalpha=0.3)
# plot!([-pi,0],-[2,2],fill=(0,:grey),fillalpha=0.3)
# plot!([2*pi,3*pi],[2,2],fill=(0,:grey),fillalpha=0.3)
# plot!([2*pi,3*pi],-[2,2],fill=(0,:grey),fillalpha=0.3)

savefig("fourier_basis.png")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier