text:
#=
Case folding for Unicode Str types

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#

_lowercase_l(ch) = ifelse(_isupper_al(ch), ch + 0x20, ch)
_uppercase_l(ch) = ifelse(_can_upper(ch),  ch - 0x20, ch)

_lowercase(ch) = is_latin(ch) ? _lowercase_l(ch) : _lowercase_u(ch)
_uppercase(ch) = is_latin(ch) ? _uppercase_l(ch) : _uppercase_u(ch)
_titlecase(ch) = is_latin(ch) ? _uppercase_l(ch) : _titlecase_u(ch)

lowercase(ch::T) where {T<:Chr} = T(_lowercase(codepoint(ch)))
uppercase(ch::T) where {T<:Chr} = T(_uppercase(codepoint(ch)))
titlecase(ch::T) where {T<:Chr} = T(_titlecase(codepoint(ch)))

lowercase(ch::ASCIIChr) = ifelse(_isupper_a(ch), ASCIIChr(ch + 0x20), ch)
uppercase(ch::ASCIIChr) = ifelse(_islower_a(ch), ASCIIChr(ch - 0x20), ch)
titlecase(ch::ASCIIChr) = uppercase(ch)

function uppercase_first(str::MaybeSub{S}) where {C<:ASCIICSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _islower_a(ch) || return str
        out = _allocate(len)
        unsafe_copyto!(out, pnt, len)
        set_codeunit!(out, ch - 0x20)
        Str(C, out)
    end
end

function lowercase_first(str::MaybeSub{S}) where {C<:ASCIICSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _isupper_a(ch) || return str
        out = _allocate(len)
        unsafe_copyto!(out, pnt, len)
        set_codeunit!(out, ch + 0x20)
        Str(C, out)
    end
end

function _upper(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:ASCIICSE}
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _islower_a(ch) && set_codeunit!(out, ch - 0x20)
        out += 1
    end
    Str(C, buf)
end

function _lower(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:ASCIICSE}
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _isupper_a(ch) && set_codeunit!(out, ch + 0x20)
        out += 1
    end
    Str(C, buf)
end

function _upper(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:LatinCSE}
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _can_upper(ch) && set_codeunit!(out, ch - 0x20)
        out += 1
    end
    Str(C, buf)
end

function uppercase(str::MaybeSub{S}) where {C<:ASCIICSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            _islower_a(get_codeunit(pnt)) && return _upper(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

function lowercase(str::MaybeSub{S}) where {C<:ASCIICSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            _isupper_a(get_codeunit(pnt)) && return _lower(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

function uppercase_first(str::MaybeSub{S}) where {C<:LatinCSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _can_upper(ch) || return str
        buf, out = _allocate(UInt8, len)
        set_codeunit!(out, ch - 0x20)
        len > 1 && unsafe_copyto!(out, pnt+1, len-1)
        Str(C, buf)
    end
end

# Special handling for characters that can't map into Latin1
function uppercase_first(str::MaybeSub{S}) where {C<:_LatinCSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        if _can_upper(ch)
            buf, out8 = _allocate(UInt8, len)
            set_codeunit!(out8, ch - 0x20)
            len > 1 && unsafe_copyto!(out8, pnt+1, len-1)
            Str(C, buf)
        elseif (ch == 0xb5) | (ch == 0xff)
            buf, out = _allocate(UInt16, len)
            set_codeunit!(out, ifelse(ch == 0xb5, 0x39c, 0x178))
            # Perform the widen operation on the rest (should be done via SIMD)
            @inbounds for i = 2:len
                set_codeunit!(out += 2, get_codeunit(pnt += 2)%UInt16)
            end
            Str(_UCS2CSE, buf)
        else
            str
        end
    end
end

function lowercase_first(str::MaybeSub{S}) where {C<:Latin_CSEs,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = pointer(str)
        ch = get_codeunit(pnt)
        _isupper(ch) || return str
        buf, out = _allocate(UInt8, len)
        set_codeunit!(out, ch + 0x20)
        len > 1 && unsafe_copyto!(out, pnt+1, len-1)
        Str(C, buf)
    end
end

lowercase(ch::T) where {T<:LatinChars} = T(_lowercase_l(codepoint(ch)))
uppercase(ch::LatinChr) = LatinChr(_uppercase_l(codepoint(ch)))

# Special handling for case where this is just an optimization of the first 256 bytes of Unicode,
# and not the 8-bit ISO 8859-1 character set
function uppercase(ch::_LatinChr)
    cb = codepoint(ch)
    _can_upper(cb) && return _LatinChr(cb - 0x20)
    # We don't uppercase 0xdf, the ß character
    cb == 0xb5 ? UCS2Chr(0x39c) : (cb == 0xff ? UCS2Chr(0x178) : ch)
end
titlecase(ch::LatinChars) = uppercase(ch)

function _upper(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:_LatinCSE}
    fin = beg + len
    cur = beg + off
    # Need to scan the rest of the string to see if _widenupper needs to be called
    while cur < fin
        ((ch = get_codeunit(cur)) == 0xb5) | (ch == 0xff) && return _widenupper(beg, off, len)
        cur += 1
    end
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _can_upper(ch) && set_codeunit!(out, ch - 0x20)
        out += 1
    end
    Str(C, buf)
end

function _widen!(dst::Ptr{T}, src::Ptr{S}, fin::Ptr{S}) where {T<:CodeUnitTypes, S<:CodeUnitTypes}
    while src < fin
        set_codeunit!(dst, get_codeunit(src)%T)
        dst += sizeof(T)
        src += sizeof(S)
    end
    nothing
end
const _narrow! = _widen!  # When this is optimized for SSE/AVX/etc. instructions, will be different

function _widenupper(beg::Ptr{UInt8}, off, len)
    buf, out = _allocate(UInt16, len)
    fin = bytoff(out, len)
    cur = beg + off
    _widen!(out, beg, cur)
    out = bytoff(out, off)
    while out < fin
        ch = get_codeunit(cur)
        if ch == 0xb5
            set_codeunit!(out, 0x39c)
        elseif ch == 0xff
            set_codeunit!(out, 0x178)
        else
            set_codeunit!(out, _can_upper(ch) ? ch - 0x20 : ch)
        end
        cur += 1
        out += 2
    end
    Str(_UCS2CSE, buf)
end

function uppercase(str::MaybeSub{S}) where {C<:LatinCSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            _can_upper(get_codeunit(pnt)) && return _upper(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

function uppercase(str::MaybeSub{S}) where {C<:_LatinCSE,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            ch = get_codeunit(pnt)
            ((ch == 0xb5) | (ch == 0xff)) && return _widenupper(beg, pnt-beg, len)
            _can_upper(ch) && return _upper(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

function _lower(::Type{C}, beg::Ptr{UInt8}, off, len) where {C<:Latin_CSEs}
    buf, out = _allocate(UInt8, len)
    fin = out + len
    unsafe_copyto!(out, beg, len)
    out += off
    while out < fin
        ch = get_codeunit(out)
        _isupper_al(ch) && set_codeunit!(out, ch + 0x20)
        out += 1
    end
    Str(C, buf)
end

function lowercase(str::MaybeSub{S}) where {C<:Latin_CSEs,S<:Str{C}}
    (len = ncodeunits(str)) == 0 && return str
    @preserve str begin
        pnt = beg = pointer(str)
        fin = beg + len
        while pnt < fin
            _isupper_al(get_codeunit(pnt)) && return _lower(C, beg, pnt-beg, len)
            pnt += 1
        end
    end
    str
end

_can_upper_latin(ch)      = _can_upper(ch) | (ch == 0xb5) | (ch == 0xff)
_can_upper_only_latin(ch) = _can_upper_l(ch) | (ch == 0xb5) | (ch == 0xff)

@inline _can_upper_ch(ch) =
    ch <= 0x7f ? _islower_a(ch) : (ch <= 0xff ? _can_upper_only_latin(ch) : _islower_u(ch))
@inline _can_lower_ch(ch) =
    ch <= 0x7f ? _isupper_a(ch) : (ch <= 0xff ? _isupper_l(ch) : _isupper_u(ch))

# result must have at least one character > 0xff, so if the only character(s)
# > 0xff became <= 0xff, then the result may need to be narrowed and returned as _LatinStr

function _lower(::Type{C}, beg, off, len) where {C<:_UCS2CSE}
    CU = codeunit(C)
    buf, out = _allocate(CU, len)
    unsafe_copyto!(out, beg, len)
    fin = out + (len*sizeof(CU))
    out += off
    flg = false
    while out < fin
        ch = get_codeunit(out)
        if ch <= 0x7f
            _isupper_a(ch) && set_codeunit!(out, ch += 0x20)
        elseif ch <= 0xff
            _isupper_l(ch) && set_codeunit!(out, ch += 0x20)
        elseif _isupper_u(ch)
            ch = _lowercase_u(ch)
            flg = ch <= 0xff
            set_codeunit!(out, ch)
        end
        out += sizeof(CU)
    end
    if flg && is_latin(buf)
        out = pointer(buf)
        buf = _allocate(len)
        _narrow!(pointer(buf), out, out + len)
        Str(_LatinCSE, buf)
    else
        Str(C, buf)
    end
end

function _lower(::Type{C}, beg, off, len) where {C<:Union{UCS2CSE,UTF32_CSEs}}
    CU = codeunit(C)
    buf, out = _allocate(CU, len)
    unsafe_copyto!(out, beg, len)
    fin = out + (len*sizeof(CU))
    out += off
    while out < fin
        ch = get_codeunit(out)
        if ch <= 0x7f
            _isupper_a(ch) && set_codeunit!(out, ch += 0x20)
        elseif ch <= 0xff
            _isupper_l(ch) && set_codeunit!(out, ch += 0x20)
        elseif _isupper_u(ch)
            set_codeunit!(out, _lowercase_u(ch))
        end
        out += sizeof(CU)
    end
    Str(C, buf)
end

function lowercase(str::MaybeSub{S}) where {C<:Union{UCS2_CSEs,UTF32_CSEs},S<:Str{C}}
    @preserve str begin
        CU = codeunit(C)
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            _can_lower_ch(get_codeunit(pnt)) && return _lower(C, beg, pnt-beg, ncodeunits(str))
            pnt += sizeof(CU)
        end
    end
    str
end

function _upper(::Type{C}, beg, off, len) where {C<:Union{UCS2_CSEs,UTF32_CSEs}}
    CU = codeunit(C)
    buf, out = _allocate(CU, len)
    unsafe_copyto!(out, beg, len)
    fin = out + (len*sizeof(CU))
    out += off
    while out < fin
        ch = get_codeunit(out)
        if ch <= 0x7f
            _islower_a(ch) && set_codeunit!(out, ch -= 0x20)
        elseif ch > 0xff
            _islower_u(ch) && set_codeunit!(out, _uppercase_u(ch))
        elseif _can_upper(ch)
            set_codeunit!(out, ch -= 0x20)
        elseif ch == 0xb5
            set_codeunit!(out, 0x39c)
        elseif ch == 0xff
            set_codeunit!(out, 0x178)
        end
        out += sizeof(CU)
    end
    Str(C, buf)
end

function uppercase(str::MaybeSub{S}) where {C<:Union{UCS2_CSEs,UTF32_CSEs},S<:Str{C}}
    @preserve str begin
        CU = codeunit(C)
        pnt = beg = pointer(str)
        fin = beg + sizeof(str)
        while pnt < fin
            _can_upper_ch(get_codeunit(pnt)) && return _upper(C, beg, pnt-beg, ncodeunits(str))
            pnt += sizeof(CU)
        end
        str
    end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct EmojiApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_emojiList(_api::EmojiApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/emoji.list", ["slackAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded"] : [_mediaType])
    return _ctx
end

"""

Lists custom emoji for a team.
Return: Dict{String, Any}
"""
function emojiList(_api::EmojiApi; _mediaType=nothing)
    _ctx = _swaggerinternal_emojiList(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function emojiList(_api::EmojiApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_emojiList(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export emojiList

ext: jl


text:
function test_forwarddiff_1()

    (Q, T, H) = (3, 2000, 40)
    p = randn(1, T)
    r = randn(H, Q)
    f1 = x -> begin
        y = p .* x
        A = r * y
        B = PERK.A_mul_A′(A)
        norm(B)
    end
    f2 = x -> begin
        y = p .* x
        A = r * y
        B = A * A'
        norm(B)
    end

    x = randn(Q)
    g1 = similar(x)
    g2 = similar(x)
    config1 = ForwardDiff.GradientConfig(f1, x, ForwardDiff.Chunk{Q}())
    config2 = ForwardDiff.GradientConfig(f2, x, ForwardDiff.Chunk{Q}())
    ForwardDiff.gradient!(g1, f1, x, config1)
    ForwardDiff.gradient!(g2, f2, x, config2)

    return g1 ≈ g2

end

@testset "ForwardDiff" begin

    @test test_forwarddiff_1()

end


text:
"""
    GutzwillerSampling(::AbstractHamiltonian; g)

Wrapper over any [`AbstractHamiltonian`](@ref) that implements Gutzwiller sampling. In this
importance sampling scheme the Hamiltonian is modified as follows.

```math
\\tilde{H}_{ij} = H_{ij} e^{-g(H_{jj} - H_{ii})}
```

# Constructor

* `GutzwillerSampling(::AbstractHamiltonian, g)`
* `GutzwillerSampling(::AbstractHamiltonian; g)`

After construction, we can access the underlying hamiltonian with `G.hamiltonian` and the
`g` parameter with `G.g`.

# Example

```jldoctest
julia> H = HubbardMom1D(BoseFS{3}((1,1,1)); u=6.0, t=1.0)
HubbardMom1D(BoseFS{3,3}((1, 1, 1)); u=6.0, t=1.0)

julia> I = GutzwillerSampling(H, g=0.3)
GutzwillerSampling(HubbardMom1D(BoseFS{3,3}((1, 1, 1)); u=6.0, t=1.0); g=0.3)

julia> get_offdiagonal(H, BoseFS((2, 1, 0)), 1)
(BoseFS{3,3}((1, 0, 2)), 2.0)

julia> get_offdiagonal(I, BoseFS((2, 1, 0)), 1)
(BoseFS{3,3}((1, 0, 2)), 0.8131393194811987)
```
"""
struct GutzwillerSampling{A,T,H<:AbstractHamiltonian{T},G} <: AbstractHamiltonian{T}
    # The A parameter sets whether this is an adjoint or not.
    # The G parameter is g parameter value.
    hamiltonian::H
end

function GutzwillerSampling(h, g)
    G = eltype(h)(g)
    return GutzwillerSampling{false,eltype(h),typeof(h),G}(h)
end
GutzwillerSampling(h; g) = GutzwillerSampling(h, g)

function Base.show(io::IO, h::GutzwillerSampling{A}) where {A}
    A && print(io, "adjoint(")
    print(io, "GutzwillerSampling(", h.hamiltonian, "; g=", h.g, ")")
    A && print(io, ")")
end

starting_address(h::GutzwillerSampling) = starting_address(h.hamiltonian)

LOStructure(::Type{<:GutzwillerSampling{<:Any,<:Any,H}}) where {H} = _lo_str(LOStructure(H))
_lo_str(::LOStructure) = AdjointKnown()
_lo_str(::AdjointUnknown) = AdjointUnknown()

function LinearAlgebra.adjoint(h::GutzwillerSampling{A,T,<:Any,G}) where {A,T,G}
    h_adj = h.hamiltonian'
    return GutzwillerSampling{!A,T,typeof(h_adj),G}(h_adj)
end

dimension(::Type{T}, h::GutzwillerSampling) where T = dimension(T, h.hamiltonian)

Base.getproperty(h::GutzwillerSampling, s::Symbol) = getproperty(h, Val(s))
Base.getproperty(h::GutzwillerSampling{<:Any,<:Any,<:Any,G}, ::Val{:g}) where G = G
Base.getproperty(h::GutzwillerSampling, ::Val{:hamiltonian}) = getfield(h, :hamiltonian)

# Forward some interface functions.
num_offdiagonals(h::GutzwillerSampling, add) = num_offdiagonals(h.hamiltonian, add)
diagonal_element(h::GutzwillerSampling, add) = diagonal_element(h.hamiltonian, add)

function gutzwiller_modify(matrix_element, is_adjoint, g, diag1, diag2)
    if is_adjoint
        return matrix_element * exp(-g * (diag1 - diag2))
    else
        return matrix_element * exp(-g * (diag2 - diag1))
    end
end

function get_offdiagonal(h::GutzwillerSampling{A}, add1, chosen) where A
    add2, matrix_element = get_offdiagonal(h.hamiltonian, add1, chosen)
    diag1 = diagonal_element(h, add1)
    diag2 = diagonal_element(h, add2)
    return add2, gutzwiller_modify(matrix_element, A, h.g, diag1, diag2)
end

struct GutzwillerOffdiagonals{
    F,T,A,G,H<:AbstractHamiltonian,N<:AbstractOffdiagonals{F,T}
}<:AbstractOffdiagonals{F,T}
    hamiltonian::H
    diag::T
    offdiagonals::N
end

function offdiagonals(ham::GutzwillerSampling{A,<:Any,<:Any,G}, a) where {A,G}
    hps = offdiagonals(ham.hamiltonian, a)
    diag = diagonal_element(ham, a)
    return GutzwillerOffdiagonals{typeof(a),eltype(ham),A,G,typeof(ham),typeof(hps)}(
        ham, diag, hps
    )
end

function Base.getindex(h::GutzwillerOffdiagonals{F,T,A,G}, i)::Tuple{F,T} where {F,T,A,G}
    add2, matrix_element = h.offdiagonals[i]
    diag2 = diagonal_element(h.hamiltonian, add2)
    return add2, gutzwiller_modify(matrix_element, A, G, h.diag, diag2)
end

Base.size(h::GutzwillerOffdiagonals) = size(h.offdiagonals)


text:
##############################################################################
##
## AbstractDataset interface
##
##############################################################################

index(ds::Dataset) = getfield(ds, :colindex)
_copy_grouping_info!(ds::Dataset, ds2::Dataset) = _copy_grouping_info!(index(ds), index(ds2))
_reset_grouping_info!(ds::Dataset) = _reset_grouping_info!(index(ds))
_attributes(ds::Dataset) = getfield(ds, :attributes)
# this function grants the access to the internal storage of columns of the
# `Dataset` and its use is unsafe. If the returned vector is mutated then
# make sure that:
# 1. `AbstractRange` columns are not added to a `Dataset`
# 2. all inserted columns use 1-based indexing
# 3. after several mutating operations on the vector are performed
#    each element (column) has the same length
# 4. if length of the vector is changed that the index of the `Dataset`
#    is adjusted appropriately

_onebased_check_error() =
    throw(ArgumentError("Currently InMemoryDatasets.jl supports only columns " *
                        "that use 1-based indexing"))
_onebased_check_error(i, col) =
    throw(ArgumentError("Currently InMemoryDatasets.jl supports only " *
                        "columns that use 1-based indexing, but " *
                        "column $i has starting index equal to $(firstindex(col))"))

# note: these type assertions are required to pass tests
nrow(ds::Dataset) = ncol(ds) > 0 ? length(_columns(ds)[1])::Int : 0
ncol(ds::Dataset) = length(index(ds))

##############################################################################
##
## Dataset consistency check
##
##############################################################################

corrupt_msg(ds::Dataset, i::Integer) =
    "Data set is corrupt: length of column " *
    ":$(_names(ds)[i]) ($(length(ds[!, i]))) " *
    "does not match length of column 1 ($(length(ds[!, 1]))). " *
    "The column vector has likely been resized unintentionally " *
    "(either directly or because it is shared with another data set)."

function _check_consistency(ds::Dataset)
    cols, idx = _columns(ds), index(ds)

    for (i, col) in enumerate(cols)
        firstindex(col) != 1 && _onebased_check_error(i, col)
    end

    ncols = length(cols)
    @assert length(idx.names) == length(idx.lookup) == ncols
    ncols == 0 && return nothing
    nrows = length(cols[1])
    for i in 2:length(cols)
        @assert length(cols[i]) == nrows corrupt_msg(ds, i)
    end
    nothing
end

function _check_consistency(ds::AbstractDataset)
    if ds isa SubDataset
        @assert length(index(ds).remap) == length(index(parent(ds))) "The parent data set which this view is based on, has been modified. To fix the issue recreate the view"
    end
    _check_consistency(parent(ds))
end

##############################################################################
##
## getindex() definitions
##
##############################################################################

# ds[SingleRowIndex, SingleColumnIndex] => Scalar
@inline function Base.getindex(ds::Dataset, row_ind::Integer,
                               col_ind::Union{Signed, Unsigned})
    cols = _columns(ds)
    @boundscheck begin
        if !checkindex(Bool, axes(cols, 1), col_ind)
            throw(BoundsError(ds, (row_ind, col_ind)))
        end
        if !checkindex(Bool, axes(ds, 1), row_ind)
            throw(BoundsError(ds, (row_ind, col_ind)))
        end
    end

    @inbounds cols[col_ind][row_ind]
end

@inline function Base.getindex(ds::Dataset, row_ind::Integer,
                               col_ind::SymbolOrString)
    selected_column = index(ds)[col_ind]
    @boundscheck if !checkindex(Bool, axes(ds, 1), row_ind)
        throw(BoundsError(ds, (row_ind, col_ind)))
    end
    @inbounds _columns(ds)[selected_column][row_ind]
end

_getindex(ds, r, c, f) = map(f, view(_columns(ds)[c], r))

# ds[MultiRowIndex, SingleColumnIndex] => AbstractVector, copy
@inline function Base.getindex(ds::Dataset, row_inds::AbstractVector, col_ind::ColumnIndex; mapformats = false)
    selected_column = index(ds)[col_ind]
    @boundscheck if !checkindex(Bool, axes(ds, 1), row_inds)
        throw(BoundsError(ds, (row_inds, col_ind)))
    end
    if mapformats
        _getindex(ds, row_inds, col_ind, getformat(ds, col_ind))
    else
        @inbounds return _columns(ds)[selected_column][row_inds]
    end
end

@inline function Base.getindex(ds::Dataset, row_inds::Not, col_ind::ColumnIndex; mapformats = false)
    if mapformats
        _getindex(ds, axes(ds, 1)[row_inds], col_ind, getformat(ds, col_ind))
    else
        ds[axes(ds, 1)[row_inds], col_ind]
    end
end

# ds[:, SingleColumnIndex] => AbstractVector
function Base.getindex(ds::Dataset, row_inds::Colon, col_ind::ColumnIndex; mapformats = false)
    selected_column = index(ds)[col_ind]
    if mapformats
        map(getformat(ds, col_ind), _columns(ds)[selected_column])
    else
        copy(_columns(ds)[selected_column])
    end
end

# ds[!, SingleColumnIndex] => AbstractVector, the same vector
@inline function Base.getindex(ds::Dataset, ::typeof(!), col_ind::Union{Signed, Unsigned})
    cols = _columns(ds)
    @boundscheck if !checkindex(Bool, axes(cols, 1), col_ind)
        throw(BoundsError(ds, (!, col_ind)))
    end
    @inbounds DatasetColumn(col_ind, ds, _columns(ds)[col_ind])
end

function Base.getindex(ds::Dataset, ::typeof(!), col_ind::SymbolOrString)
    selected_column = index(ds)[col_ind]
    ds[!, selected_column]
end

# ds[MultiRowIndex, MultiColumnIndex] => Dataset
function _threaded_permute(x, perm; threads = true)
    x_cpy = similar(x, 1)
    resize!(x_cpy, length(perm))
    @_threadsfor threads for i in 1:length(x_cpy)
        x_cpy[i] = x[perm[i]]
    end
    x_cpy
end

# Create Dataset
function _threaded_getindex(selected_rows::AbstractVector,
    selected_columns::AbstractVector,
    ds_columns::AbstractVector,
    idx::AbstractIndex)
    # FIXME threading should be done along rows rather than columns
    # @static if VERSION >= v"1.4"
    #     if length(selected_rows) >= 1_000_000 && Threads.nthreads() > 1
    #         new_columns = Vector{AbstractVector}(undef, length(selected_columns))
    #         @sync for i in eachindex(new_columns)
    #             Threads.@spawn new_columns[i] = ds_columns[selected_columns[i]][selected_rows]
    #         end
    #         return Dataset(new_columns, idx, copycols=false)
    #     else
    #         return Dataset(AbstractVector[ds_columns[i][selected_rows] for i in selected_columns],
    #                          idx, copycols=false)
    #     end
    # else
    #     return Dataset(AbstractVector[ds_columns[i][selected_rows] for i in selected_columns],
    #                      idx, copycols=false)
    # end
    new_columns = Vector{AbstractVector}(undef, length(selected_columns))
    # for many columns threads over columns
    if length(selected_columns) > 100
        Threads.@threads for j in 1:length(selected_columns)
            new_columns[j] = ds_columns[selected_columns[j]][selected_rows]
        end
    else
        for j in 1:length(selected_columns)
            # TODO should we be careful about types string, number, PooledArrays???
            if DataAPI.refpool(ds_columns[selected_columns[j]]) !== nothing
                pa = ds_columns[selected_columns[j]]
                if pa isa PooledArray
                    # we could use copy but it will be inefficient for small selected_rows
                    new_columns[j] = PooledArray(PooledArrays.RefArray(_threaded_permute(pa.refs, selected_rows)), DataAPI.invrefpool(pa), DataAPI.refpool(pa), PooledArrays.refcount(pa))
                else
                    # for other pooled data(like Categorical arrays) we don't have optimised path
                    new_columns[j] = pa[selected_rows]
                end
            else
                new_columns[j] = _threaded_permute(ds_columns[selected_columns[j]], selected_rows)
            end
        end
    end
    return Dataset(new_columns, idx, copycols=false)
end

# Create Dataset
@inline function Base.getindex(ds::Dataset, row_inds::AbstractVector{T},
                               col_inds::MultiColumnIndex) where T
    @boundscheck if !checkindex(Bool, axes(ds, 1), row_inds)
        throw(BoundsError(ds, (row_inds, col_inds)))
    end
    selected_columns = index(ds)[col_inds]

    u = _names(ds)[selected_columns]
    lookup = Dict{Symbol, Int}(zip(u, 1:length(u)))
    dsformat = index(ds).format
    format = Dict{Int, Function}()
    for i in 1:length(selected_columns)
      if haskey(dsformat, selected_columns[i])
        push!(format, i => dsformat[selected_columns[i]])
      end
    end
    # use this constructor to avoid checking twice if column names are not
    # duplicate as index(ds)[col_inds] already checks this
    idx = Index(lookup, u, format)

    if length(selected_columns) == 1
        newds = Dataset(AbstractVector[_columns(ds)[selected_columns[1]][row_inds]],
                         idx, copycols=false)
        setinfo!(newds, _attributes(ds).meta.info[])
        newds
    else
        # Computing integer indices once for all columns is faster
        selected_rows = T === Bool ? _findall(row_inds) : row_inds
        newds = _threaded_getindex(selected_rows, selected_columns, _columns(ds), idx)
        setinfo!(newds, _attributes(ds).meta.info[])
        newds
    end
end

# Create Dataset
@inline function Base.getindex(ds::Dataset, row_inds::AbstractVector{T}, ::Colon) where T
    @boundscheck if !checkindex(Bool, axes(ds, 1), row_inds)
        throw(BoundsError(ds, (row_inds, :)))
    end
    dsidx = index(ds)
    idx = Index(copy(dsidx.lookup), copy(dsidx.names), copy(dsidx.format))

    if ncol(ds) == 1
        newds = Dataset(AbstractVector[_columns(ds)[1][row_inds]], idx, copycols=false)
        setinfo!(newds, _attributes(ds).meta.info[])
        newds
    else
        # Computing integer indices once for all columns is faster
        selected_rows = T === Bool ? _findall(row_inds) : row_inds
        newds = _threaded_getindex(selected_rows, 1:ncol(ds), _columns(ds), idx)
        setinfo!(newds, _attributes(ds).meta.info[])
        newds
    end
end

# Create Dataset
@inline Base.getindex(ds::Dataset, row_inds::Not, col_inds::MultiColumnIndex) =
    ds[axes(ds, 1)[row_inds], col_inds]

# ds[:, MultiColumnIndex] => Dataset
# Create Dataset
# In select we should check if the attributes can be passed and also check how to set info
Base.getindex(ds::Dataset, row_ind::Colon, col_inds::MultiColumnIndex) =
    select(ds, col_inds)

# ds[!, MultiColumnIndex] => Dataset
# Create Dataset
# the same as : we should check about attributes which can be passed
function Base.getindex(ds::Dataset, row_ind::typeof(!), col_inds::MultiColumnIndex)
    colsidx = index(ds)[col_inds]
    view(ds, :, colsidx)
end
    # throw(ArgumentError("syntax `ds[!, col]` is only valid for single column, for multiple columns use `ds[:, cols]` or `select(ds, col...)`"))


text:
using BigNums: BITS, normalize!, set_one!, set_zero!, bitneg!

### Properties ###

###################
### Operational ###
###################

begin
    function checkOne()
        a = one(ArbUInt)
        a.data == [one(ArbDigit)]
    end

    function checkZero()
        a = zero(ArbUInt)
        isempty(a.data)
    end

    function isnormalized(a::ArbUInt)
        findlast(iszero, a.data) != lastindex(a.data)
    end

    function constructNormalized(v::Vector{ArbDigit})
        isnormalized(ArbUInt(v))
    end

    function constructNormalized(v::Union{ArbDigit, DoubleArbDigit})
        isnormalized(ArbUInt(v))
    end

    function setOne(a)
        set_one!(deepcopy(a)) == one(a)
    end

    function setZero(a)
        set_zero!(deepcopy(a)) == zero(a)
    end
end

################
### Ordering ###
################

begin
    function orderTransitive(a,b,c)
        if a < b && b < c
            return a < c
        elseif a < c && c < b
            return a < b
        elseif b < c && c < a
            return b < a
        elseif b < a && a < c
            return b < c
        elseif c < a && a < b
            return c < b
        elseif c < b && b < a
            return c < a
        else
            return a == b || b == c || a == c
        end
    end

    function orderReversal(a,b)
        if a < b
            return b > a
        elseif b < a
            return a > b
        else
            return a == b
        end
    end

    function trichotomy(a,b)
        x1 = a < b
        x2 = a == b
        x3 = a > b
        if x1
            return x2 == x3 == false
        elseif x2
            return x1 == x3 == false
        elseif x3
            return x1 == x2 == false
        end
    end

    function orderPreservedAddition((small,big),constant)
        (small + constant) < (big + constant)
    end
end

##################
### Arithmetic ###
##################

begin ## Addition ##
    function commutative_add(a::ArbUInt, b::ArbUInt)
        c = a + b
        d = b + a
        c == d
    end

    function associative_add(a::ArbUInt, b::ArbUInt, c::ArbUInt)
        d = a + (b + c)
        e = (a + b) + c
        d == e
    end

    function identity_add(a::ArbUInt)
        a == (a + zero(ArbUInt)) && a == (zero(ArbUInt) + a)
    end
end

begin ## Multiplication ##
    function long_karatsuba_eq(a::AbstractVector{ArbDigit}, b::AbstractVector{ArbDigit})
        len = length(a) + length(b)
        long_buf = zeros(ArbDigit, len)
        kara_buf = zeros(ArbDigit, len)

        short, long = length(a) < length(b) ? (a,b) : (b,a)

        BigNums.long_mul!(long_buf, long, short)
        BigNums.karatsuba!(kara_buf, long, short)

        all(Base.splat(==), zip(kara_buf, long_buf))
    end

    function commutative_mul(a::ArbUInt, b::ArbUInt)
        c = a * b
        d = b * a
        c == d
    end

    function associative_mul(a::ArbUInt, b::ArbUInt, c::ArbUInt)
        d = a * (b * c)
        e = (a * b) * c
        d == e
    end

    function distributive_mul(a::ArbUInt, b::ArbUInt, c::ArbUInt)
        d = a * (b + c)
        e = (a*b) + (a*c)
        d == e
    end

    function identity_mul(a::ArbUInt)
        a == a * one(ArbUInt) && a == one(ArbUInt) * a
    end

    function zero_mul(a::ArbUInt)
        iszero(a * zero(ArbUInt)) && iszero(zero(ArbUInt) * a)
    end
end

###############
### Bitwise ###
###############

commutative(op) = (a,b) -> commutative(a,b,op)
function commutative(a, b, op)
    op(a,b) == op(b,a)
end

begin ## Shifts ##
    function leftShiftTrailLeadCorrect(a, sh)
        b = a << sh
        iszero(a) && return trailing_zeros(a) == trailing_zeros(b) == 0

        trail_correct = trailing_zeros(a) + sh == trailing_zeros(b)

        sh %= BITS # => sh < BITS (the part we care about)

        lead_a = leading_zeros(a)
        lead_b = leading_zeros(b)
        lead_correct = begin
            if iszero(sh) # did we just shift the top most word up a number of digits?
                lead_a == lead_b
            elseif sh > lead_a # ok, only a part word overall, but did we overflow?
                lead_b == BITS - (sh - lead_a)
            elseif sh < lead_a # what about without overflow?
                lead_b == lead_a - sh
            else # and what if we shift exactly up to the boundary?
                iszero(lead_b)
            end
        end

        return trail_correct && lead_correct
    end

    function leftShiftNoDataLoss(a, sh)
        b = a << sh
        count_ones(b) == count_ones(a)
    end

    function additionOfShifts(a, (sh1, sh2)::Tuple{T,T}) where T <: Unsigned
        T1 = widen(T)
        ((a << sh1) << sh2) == (a << (T1(sh1) + T1(sh2)))
    end

    function shiftAntiInverse(a, sh)
        ((a << sh) >> sh) == a
    end

    function shiftValueSubtraction(a, (sh1,sh2)::Tuple{T,T}) where T <: Unsigned
        # property assumes sh1 <= sh2
        b = (a << sh2) >> sh1
        c = a << (sh2 - sh1)
        b == c
    end

    function rightShiftSame(a::ArbUInt, sh)
        (a >>> sh) == (a >> sh)
    end

    function shiftIdentity(a)
        (a >> 0x0) == a && (a << 0x0) == a
    end
end

begin ## Negation ##
    function inverse(a)
        a == ~(~a)
    end
end

begin ## And ##
    function zeroing_and(a)
        (a & zero(a)) == zero(a) && (zero(a) & a) == zero(a)
    end

    function identity_and(a)
        (a & a) == a
    end
end

begin ## Or ##
    function identity_or(a, b)
       (a | a) == a
    end
end

begin ## Xor ##
    function zeroing_xor(a)
        (a ⊻ a) == zero(a)
    end

    function identity_xor(a)
        (a ⊻ zero(a)) == a
    end
end

begin ## Nor ##
    function chaining_nor(a,b)
        nor(a,b) == normalize!(~(a | b))
    end
end

begin ## Nand ##
    function chaining_nand(a,b)
        nand(a,b) == normalize!(~(a & b))
    end
end

begin ## Combined Properties ##
    function lossyShift(a, sh)
        b = (a >> sh) << sh
        # mask out the parts we shifted out in b
        # FIXME: This generation of the mask is broken for when we have more than one word
        mask = ArbUInt(unsigned(-1) % ArbDigit) << sh
        # mask = (one(ArbUInt) << sh) - 0x1 # this is proper, but we don't have subtraction yet
        pop!(mask.data)
        bitneg!(mask)
        (b ⊻ (a & mask)) == a
    end
end

begin #= Stringification & Parsing =#
    function parsestring(str, base)
        string(parse(ArbUInt, str); base=base, pad=length(str)) == str
    end

    function stringparse(a, base)
        parse(ArbUInt, string(a; base=base)) == a
    end

    function throwNonHex(str)
        try
            parse(ArbUInt, str; base=16)
        catch ex
            ex isa ArgumentError && return true
            rethrow(ex)
        end
        return false
    end
end


text:
function zad35(d, el1, el2)
    A = zeros(d,d)
    A[:, 1:2:end] .= el1
    A[:, 2:2:end] .= el2
    return A

end


text:
export convCircKernel ,getConvCircKernel
## For the functions nImgIn, nImgOut, nFeatIn, nFeatOut, nTheta, getOp, initTheta : see AbstractConvKernel.jl
## All convKernel types are assumed to have fields nImage and sK
mutable struct convCircKernel{T} <: AbstractConvKernel{T}
    nImg :: Array{Int,1}
    sK   :: Array{Int,1}
    iK   :: Array{Int,1}
end

function nImgOut(this::convCircKernel)
    return [this.nImg[1]; this.nImg[2]; this.sK[3]]
 end

function getConvCircKernel(TYPE::Type,nImg,sK)

    n3 = [nImg[1],nImg[2],sK[3]]
    Ind = zeros(n3[1],n3[2],n3[3])
    mid = div.(n3,2)
    midK = div.(sK[1:2],2)
    Ind[mid[1]-midK[1]:mid[1]+midK[1],
        mid[2]-midK[2]:mid[2]+midK[2],:] = reshape(collect(1:prod(sK)),sK[1],sK[2],sK[3])
    Ind = circshift(Ind,1-mid)

    iK = find(Ind)
	return convCircKernel{TYPE}(nImg,sK,iK)
end

export Amv
function Amv(this::convCircKernel{T},theta::Array{T},Yin::Array{T}) where {T<:Number}

    nex   = div(length(Y),prod(nImgIn(this)))

    n     = nImgIn(this);

    Y     = reshape(Yin,(n[1],n[2],this.sK[3],nex))
    # compute convolution
    K = zeros(n[1],n[2],this.sK[3])

    K[this.iK] = theta[:];

    Kh = fft(K,(1,2,3))
    Yh = fft(Y,(1,2,3))

    Z  = real(ifft(Kh.*Yh,(3,2,1)))
    return Z
end

export ATmv
function ATmv(this::convCircKernel{T},theta::Array{T},Zin::Array{T}) where {T<:Number}

    nex   =  div(numel(Z),prod(nImgIn(this)))
    n     = nImgIn(this)
    Z     = reshape(Zin,(n[1],n[2],this.sK[3],nex));
    # compute convolution T
    K = zeros(n[1],n[2],this.sK[3])
    K[this.iK] = theta

    Kh = fft(K,(1,2,3))
    Zh = fft(Z,(1,2,3))

    Y  = real(ifft(conj(Kh).*Zh,(3,2,1)))
    return Y
end

function Jthetamv(this::convCircKernel{T},dtheta::Array{T},dummy::Array{T},Yin::Array{T},temp=nothing) where {T<:Number}

    nex    =  div(numel(Y),nFeatIn(this))
    Y      = reshape(Yin,:,nex)
    Z      = Amv(this,dtheta,Y)
    return Z
end

function JthetaTmv(this::convCircKernel{T},Zin::Array{T},dummy::Array{T},Yin::Array{T}) where {T<:Number}
    #  derivative of Z*(A(theta)*Y) w.r.t. theta

    nex  =  div(numel(Y),nFeatIn(this))

    n  = nImgIn(this)

    Z     = reshape(Zin,n[1],n[2],this.sK[3],nex)
    Y     = reshape(Yin,n[1],n[2],this.sK[3],nex)

    Zh = fft(Z,(1,2,3))
    Yh = fft(Y,(1,2,3))
    YZ = sum(conj(Yh).*Zh,4)
    K  = real(ifft(YZ,(3,2,1)))
    dtheta = K[this.iK];

end


text:
push!(LOAD_PATH, "/Users/eroesch/Documents/phd/neural-ode/")
using Flux, DiffEqFlux, DifferentialEquations, Plots

## Setup ODE to optimize
function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(lotka_volterra,u0,tspan,p)

# Verify ODE solution
sol = solve(prob,Tsit5())
plot(sol)

# Generate data from the ODE
sol = solve(prob,Tsit5(),saveat=0.1)
A = sol[1,:] # length 101 vector
t = 0:0.1:10.0
scatter!(t,A)

# Build a neural network that sets the cost as the difference from the
# generated data and 1

p = param([2.2, 1.0, 2.0, 0.4]) # Initial Parameter Vector
function predict_rd() # Our 1-layer neural network
  diffeq_rd(p,prob,Tsit5(),saveat=0.1)[1,:]
end
loss_rd() = sum(abs2,x-1 for x in predict_rd()) # loss function

# Optimize the parameters so the ODE's solution stays near 1

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  display(loss_rd())
  # using `remake` to re-create our `prob` with current parameters `p`
  display(plot(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6)))
end
# Display the ODE with the initial parameter values.
cb()
Flux.train!(loss_rd, [p], data, opt, cb = cb)

max_stars_repo_name: LislPisl/model-zoo

max_issues_repo_name: LislPisl/model-zoo

max_forks_repo_name: LislPisl/model-zoo


text:
module Ex7_3
using DescriptorSystems
using LinearAlgebra
using Polynomials
using Test

# Example 7.3 - Nullspace-based synthesis
# Uses only DescriptorSystems.jl
println("Example 7.3")

A = diagm([ -2, -3, -1, 2, 3]);
Bu = [1 1 0 0 0]'; Bd = [0 0 2 0 0]';
Bf = [0 0 0 2 2;0 0 0 0 0]';
C = [-1 0 -1 1.5 0; 0 -1 0 0 2.5];
Du = [ 1 1]'; Dd = [1 0]'; Df = [1 0; 1 1];
p = 2; mu = 1; md = 1; mf = 2;       # enter dimensions
sys = dss(A,[Bu Bd Bf],C,[Du Dd Df]); # define system

# compute [Q Rf], where Q is a left nullspace basis of
# [Gu Gd;I 0] and Rf = Q*[Gf;0];
Q_Rf = glnull([sys;eye(mu,mu+md+mf)],mf)[1];

# stabilize using left coprime factorization
Q_Rf, M = glcf(Q_Rf, sdeg = -3, smarg = -2);

# compute Q and Rf; scale to meet example
Q = -sqrt(2)*Q_Rf[:,1:p+mu]; Rf = -sqrt(2)*Q_Rf[:,p+mu+1:end];

# display results
println("Q = $(dss2rm(Q))")
println("Rf = $(dss2rm(Rf))")
println("M = $(dss2rm(M))")

end # module

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
using RegistryCLI: main
using LibGit2
using Test

function grep(str, file)
    contents = read(file) |> String
    findfirst(str, contents) !== nothing
end

function count_commits(path)
    repo   = LibGit2.GitRepo(path)
    walker = LibGit2.GitRevWalker(repo)
    LibGit2.count(walker) do _...
        true
    end
end

@testset "RegistryCLI" begin
    main("-h")

    Base.Filesystem.mktempdir() do tmpdir
        reg_url  = joinpath(tmpdir, "registry_public.git")  # Public registry URL
        reg_path = joinpath(tmpdir, "registry_local")       # Path to local clone

        # Create git repo for the public Registry
        LibGit2.init(#=path=# reg_url,
                     #=bare=# true)

        # Initialize public Registry
        main("create", "--help")
        main("-v", "create", reg_url, reg_path,
             "--description", "Test registry")

        @test isfile(joinpath(reg_path, "Registry.toml"))
        @test count_commits(reg_path) == 1

        # Register package
        main("add", "--help")
        cd(reg_path) do
            main("-v", "add", "https://github.com/JuliaLang/Example.jl.git",
                 "--ref", "v0.5.3")
        end

        versions = joinpath(reg_path, "E", "Example", "Versions.toml")
        @test isfile(versions)
        @test grep("""["0.5.3"]""", versions)
        @test count_commits(reg_path) == 2

        cd(reg_path) do
            main("-v", "add", "https://github.com/JuliaLang/Example.jl.git",
                 "--ref", "29aa1b4") # v0.5.3, again
        end

        @test count_commits(reg_path) == 2

        cd(reg_path) do
            main("-v", "add", "https://github.com/JuliaLang/Example.jl.git")
        end

        @test count_commits(reg_path) == 3

        cd(reg_path) do
            main("-v", "add", "https://github.com/JuliaLang/Example.jl.git",
                 "--ref", "master") # latest version, again
        end

        @test count_commits(reg_path) == 3
    end
end


text:
module Callbacks

using CUDAapi
using Dates
using KernelAbstractions
using LinearAlgebra
using MPI
using Printf
using Statistics

using CLIMAParameters
using CLIMAParameters.Planet: day

using ..Courant
using ..Checkpoint
using ..DGMethods: courant, vars_state_conservative, vars_state_auxiliary
using ..Diagnostics
using ..GenericCallbacks
using ..MPIStateArrays
using ..ODESolvers
using ..TicToc
using ..VariableTemplates
using ..VTK
using ..Mesh.Grids: HorizontalDirection, VerticalDirection

using CuArrays, CuArrays.CUDAdrv, CuArrays.CUDAnative

_sync_device(::Type{CuArray}) = synchronize()
_sync_device(::Type{Array}) = nothing

"""
    show_updates(show_updates_opt, solver_config, user_info_callback)

Return a callback function that shows simulation updates at the specified
interval and also invokes the user-specified info callback.
"""
function show_updates(show_updates_opt, solver_config, user_info_callback)
    timeend = solver_config.timeend

    cb_constr = CB_constructor(show_updates_opt, solver_config)
    if cb_constr !== nothing
        # set up the information callback
        upd_starttime = Ref(now())
        cb_updates = cb_constr() do (init = false)
            if init
                upd_starttime[] = now()
            else
                runtime = Dates.format(
                    convert(Dates.DateTime, Dates.now() - upd_starttime[]),
                    Dates.dateformat"HH:MM:SS",
                )
                energy = norm(solver_config.Q)
                @info @sprintf(
                    """
Update
    simtime = %8.2f / %8.2f
    runtime = %s
    norm(Q) = %.16e""",
                    ODESolvers.gettime(solver_config.solver),
                    solver_config.timeend,
                    runtime,
                    energy,
                )
                isnan(energy) && error("norm(Q) is NaN")
            end
            user_info_callback(init)
        end
        return cb_updates
    else
        return nothing
    end
end

"""
    diagnostics(diagnostics_opt, solver_config, dgn_starttime, diagnostics_config)

Return callback functions that shows simulation updates at the specified
interval and also invokes the user-specified info callback.
"""
function diagnostics(
    diagnostics_opt,
    solver_config,
    dgn_starttime,
    diagnostics_config,
)
    if diagnostics_opt != "never" && diagnostics_config !== nothing
        # set up a callback for each diagnostics group
        dgncbs = ()
        for dgngrp in diagnostics_config.groups
            cb_constr =
                CB_constructor(diagnostics_opt, solver_config, dgngrp.interval)
            cb_constr === nothing && continue
            fn = cb_constr() do (init = false)
                @tic diagnostics
                currtime = ODESolvers.gettime(solver_config.solver)
                @info @sprintf(
                    """
Diagnostics: %s
    %s at %s""",
                    dgngrp.name,
                    init ? "initializing" : "collecting",
                    string(currtime),
                )
                dgngrp(currtime, init = init)
                @toc diagnostics
                nothing
            end
            dgncbs = (dgncbs..., fn)
        end
        return dgncbs
    else
        return nothing
    end
end

"""
    vtk(vtk_opt, solver_config)

Return a callback that saves state and auxiliary variables to a VTK
file.
"""
function vtk(vtk_opt, solver_config, output_dir)
    cb_constr = CB_constructor(vtk_opt, solver_config)
    if cb_constr !== nothing
        vtknum = Ref(1)

        mpicomm = solver_config.mpicomm
        dg = solver_config.dg
        bl = dg.balance_law
        Q = solver_config.Q
        FT = eltype(Q)

        cb_vtk = cb_constr() do (init = false)
            @tic vtk
            vprefix = @sprintf(
                "%s_mpirank%04d_num%04d",
                solver_config.name,
                MPI.Comm_rank(mpicomm),
                vtknum[],
            )
            outprefix = joinpath(output_dir, vprefix)

            statenames = flattenednames(vars_state_conservative(bl, FT))
            auxnames = flattenednames(vars_state_auxiliary(bl, FT))

            writevtk(outprefix, Q, dg, statenames, dg.state_auxiliary, auxnames)

            # Generate the pvtu file for these vtk files
            if MPI.Comm_rank(mpicomm) == 0
                # name of the pvtu file
                pprefix = @sprintf("%s_num%04d", solver_config.name, vtknum[])
                pvtuprefix = joinpath(output_dir, pprefix)

                # name of each of the ranks vtk files
                prefixes = ntuple(MPI.Comm_size(mpicomm)) do i
                    @sprintf(
                        "%s_mpirank%04d_num%04d",
                        solver_config.name,
                        i - 1,
                        vtknum[],
                    )
                end
                writepvtu(pvtuprefix, prefixes, (statenames..., auxnames...))
            end

            vtknum[] += 1
            @toc vtk
            nothing
        end
        return cb_vtk
    else
        return nothing
    end
end

"""
    monitor_timestep_duration(mtd_opt, array_type, comm)

Returns a callback function that displays wall-clock time per time-step
statistics across MPI ranks in the communicator `comm`.  The times are
averaged over the `mtd_opt` time steps requested for output.  The
`array_type` is used to synchronize the compute device.
"""
function monitor_timestep_duration(mtd_opt, array_type, comm)
    if mtd_opt == "never"
        return nothing
    elseif !endswith(mtd_opt, "steps")
        @warn @sprintf(
            """
monitor-timestep-duration must be in 'steps'; %s unrecognized; disabling""",
            mtd_opt,
        )
        return nothing
    end
    steps = parse(Int, mtd_opt[1:(end - 5)])

    _sync_device(array_type)
    before = time_ns()

    cb_mtd = GenericCallbacks.EveryXSimulationSteps(steps) do
        _sync_device(array_type)
        after = time_ns()

        time_per_timesteps = after - before

        times = MPI.Gather(time_per_timesteps, 0, comm)
        if MPI.Comm_rank(comm) == 0
            ns_per_s = 1e9
            times = times ./ ns_per_s ./ steps

            @info @sprintf(
                """Wall-clock time per time-step (statistics across MPI ranks)
                   maximum (s) = %25.16e
                   minimum (s) = %25.16e
                   median  (s) = %25.16e
                   std     (s) = %25.16e
                """,
                maximum(times),
                minimum(times),
                median(times),
                std(times),
            )
        end

        _sync_device(array_type)
        before = time_ns()

        nothing
    end

    return cb_mtd
end

"""
    monitor_courant_numbers(mcn_opt, solver_config)

Return a callback function that displays Courant numbers for the simulation
at `mcn_opt` intervals.
"""
function monitor_courant_numbers(mcn_opt, solver_config)
    cb_constr = CB_constructor(mcn_opt, solver_config)
    if cb_constr !== nothing
        cb_cfl = cb_constr() do (init = false)
            Δt = solver_config.dt
            c_v = courant(
                nondiffusive_courant,
                solver_config,
                direction = VerticalDirection(),
            )
            c_h = courant(
                nondiffusive_courant,
                solver_config,
                direction = HorizontalDirection(),
            )
            ca_v = courant(
                advective_courant,
                solver_config,
                direction = VerticalDirection(),
            )
            ca_h = courant(
                advective_courant,
                solver_config,
                direction = HorizontalDirection(),
            )
            cd_v = courant(
                diffusive_courant,
                solver_config,
                direction = VerticalDirection(),
            )
            cd_h = courant(
                diffusive_courant,
                solver_config,
                direction = HorizontalDirection(),
            )
            simtime = ODESolvers.gettime(solver_config.solver)
            @info @sprintf(
                """
Courant numbers at simtime: %8.2f, Δt = %8.2f s
    Acoustic (vertical) Courant number    = %.2g
    Acoustic (horizontal) Courant number  = %.2g
    Advection (vertical) Courant number   = %.2g
    Advection (horizontal) Courant number = %.2g
    Diffusion (vertical) Courant number   = %.2g
    Diffusion (horizontal) Courant number = %.2g""",
                simtime,
                Δt,
                c_v,
                c_h,
                ca_v,
                ca_h,
                cd_v,
                cd_h,
            )
            nothing
        end
        return cb_cfl
    else
        return nothing
    end
end

"""
    checkpoint(
        checkpoint_opt,
        checkpoint_keep_one,
        solver_config,
        checkpoint_dir,
    )

Return a callback function that runs at `checkpoint_opt` intervals
and stores a simulation checkpoint into `checkpoint_dir` identified
by a running number.
"""
function checkpoint(
    checkpoint_opt,
    checkpoint_keep_one,
    solver_config,
    checkpoint_dir,
)
    cb_constr = CB_constructor(checkpoint_opt, solver_config)
    if cb_constr !== nothing
        cpnum = Ref(1)
        cb_checkpoint = cb_constr() do (init = false)
            write_checkpoint(
                solver_config,
                checkpoint_dir,
                solver_config.name,
                solver_config.mpicomm,
                cpnum[],
            )
            if checkpoint_keep_one
                rm_checkpoint(
                    checkpoint_dir,
                    solver_config.name,
                    solver_config.mpicomm,
                    cpnum[] - 1,
                )
            end
            cpnum[] += 1
            nothing
        end
        return cb_checkpoint
    else
        return nothing
    end
end

"""
    CB_constructor(interval, solver_config, default)

Parse the specified `interval` and return the appropriate `GenericCallbacks`
constructor. If `interval` is "default", then use `default` as the interval.
"""
function CB_constructor(interval::String, solver_config, default = nothing)
    mpicomm = solver_config.mpicomm
    solver = solver_config.solver
    dg = solver_config.dg
    bl = dg.balance_law
    secs_per_day = day(bl.param_set)

    if endswith(interval, "smonths")
        ticks = 30.0 * secs_per_day * parse(Float64, interval[1:(end - 7)])
        return (func) ->
            GenericCallbacks.EveryXSimulationTime(func, ticks, solver)
    elseif endswith(interval, "sdays")
        ticks = secs_per_day * parse(Float64, interval[1:(end - 5)])
        return (func) ->
            GenericCallbacks.EveryXSimulationTime(func, ticks, solver)
    elseif endswith(interval, "shours")
        ticks = 60.0 * 60.0 * parse(Float64, interval[1:(end - 6)])
        return (func) ->
            GenericCallbacks.EveryXSimulationTime(func, ticks, solver)
    elseif endswith(interval, "smins")
        ticks = 60.0 * parse(Float64, interval[1:(end - 5)])
        return (func) ->
            GenericCallbacks.EveryXSimulationTime(func, ticks, solver)
    elseif endswith(interval, "ssecs")
        ticks = parse(Float64, interval[1:(end - 5)])
        return (func) ->
            GenericCallbacks.EveryXSimulationTime(func, ticks, solver)
    elseif endswith(interval, "steps")
        steps = parse(Int, interval[1:(end - 5)])
        return (func) -> GenericCallbacks.EveryXSimulationSteps(func, steps)
    elseif endswith(interval, "hours")
        secs = 60 * 60 * parse(Int, interval[1:(end - 5)])
        return (func) ->
            GenericCallbacks.EveryXWallTimeSeconds(func, secs, mpicomm)
    elseif endswith(interval, "mins")
        secs = 60 * parse(Int, interval[1:(end - 4)])
        return (func) ->
            GenericCallbacks.EveryXWallTimeSeconds(func, secs, mpicomm)
    elseif endswith(interval, "secs")
        secs = parse(Int, interval[1:(end - 4)])
        return (func) ->
            GenericCallbacks.EveryXWallTimeSeconds(func, secs, mpicomm)
    elseif interval == "default"
        if default === nothing
            @warn "no default available; ignoring"
            return nothing
        elseif default === ""
            return nothing
        else
            return CB_constructor(default, solver_config, "")
        end
    elseif interval == "never"
        return nothing
    else
        @warn @sprintf(
            """
%s: unrecognized interval; ignoring""",
            interval,
        )
        return nothing
    end
end

function __init__()
    tictoc()
end

end # module

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{typeof(DocStringExtensions.__init__)})
    precompile(Tuple{typeof(DocStringExtensions.template_hook), LineNumberNode, Module, String, Expr})
end


text:
"""
## shape_der!

This subroutine produces derivatives of shape functions with respect
to local coordinates.

### Function
```julia
shape_der!(der, point, i)
```

### Arguments
```julia
* der::Matrix{Float64}       : Function derivative (Updated)
* points::Matrix{Float64}    : Local coordinates of integration points
* i::Int                   : Integration point
```
"""
function shape_der!(der::Matrix{Float64}, points::Matrix{Float64}, i::Int)
  #
  #   This subroutine produces derivatives of shape functions with respect
  #   to local coordinates.
  #

  ndim = size(der, 1)
  nod = size(der, 2)

  if ndim == 1                                          # 1.0 dimensional fin_els
   xi = points[i, 1]
   if nod == 2
     der[1,1] = -0.5
     der[1,2] =  0.5
   elseif nod == 3
     t1=-1.0-xi
     t2=-xi
     t3=1.0-xi
     der[1,1] = -(t3+t2)/2.0
     der[1,2] = (t3+t1)
     der[1,3] = -(t2+t1)/2.0
   elseif nod == 4
     t1=-1.0-xi
     t2=-1.0/3.0-xi
     t3=1.0/3.0-xi
     t4=1.0-xi
     der[1,1] = -(t3*t4+t2*t4+t2*t3)*9.0/16.0
     der[1,2] = (t3*t4+t1*t4+t1*t3)*27.0/16.0
     der[1,3] = -(t2*t4+t1*t4+t1*t2)*27.0/16.0
     der[1,4] = (t2*t3+t1*t3+t1*t2)*9.0/16.0
   elseif nod == 5
     t1=-1.0-xi
     t2=-0.5-xi
     t3=-xi
     t4=0.5-xi
     t5=1.0-xi
     der[1,1] = -(t3*t4*t5+t2*t4*t5+t2*t3*t5+t2*t3*t4)*2.0/3.0
     der[1,2] = (t3*t4*t5+t1*t4*t5+t1*t3*t5+t1*t3*t4)*8.0/3.0
     der[1,3] = -(t2*t4*t5+t1*t4*t5+t1*t2*t5+t1*t2*t4)*4.0
     der[1,4] = (t2*t3*t5+t1*t3*t5+t1*t2*t5+t1*t2*t3)*8.0/3.0
     der[1,5] = -(t2*t3*t4+t1*t3*t4+t1*t2*t4+t1*t2*t3)*2.0/3.0
   else
     println("Wrong number of nodes in shape_der.")
   end
  elseif ndim == 2                                      # 2.0 dimensional fin_els
   xi=points[i, 1]
   eta=points[i, 2]
   c1=xi
   c2=eta
   c3=1.0-c1-c2
   etam=0.25*(1.0-eta)
   etap=0.25*(1.0+eta)
   xim= 0.25*(1.0-xi)
   xip= 0.25*(1.0+xi)
   x2p1=2.0*xi+1.0
   x2m1=2.0*xi-1.0
   e2p1=2.0*eta+1.0
   e2m1=2.0*eta-1.0
   if nod == 3
     der[1,1] = 1.0
     der[1,3] = 0.0
     der[1,2] = -1.0
     der[2,1] = 0.0
     der[2,3] = 1.0
     der[2,2] = -1.0
   elseif nod == 6
     der[1,1] = 4.0*c1-1.0
     der[1,6] = 4.0*c2
     der[1,5] = 0.0
     der[1,4] = -4.0*c2
     der[1,3] = -(4.0*c3-1.0)
     der[1,2] = 4.0*(c3-c1)
     der[2,1] = 0.0
     der[2,6] = 4.0*c1
     der[2,5] = 4.0*c2-1.0
     der[2,4] = 4.0*(c3-c2)
     der[2,3] = -(4.0*c3-1.0)
     der[2,2] = -4.0*c1
   elseif nod == 10
     der[1,1] = (27.0*c1*c1-18.0*c1+2.0)/2.0
     der[1,9] = (9.0*(6.0*c1-1.0)*c2)/2.0
     der[1,8] = (9.0*(3.0*c2-1.0)*c2)/2.0
     der[1,7] = 0.0
     der[1,6] = -(9.0*(3.0*c2-1.0)*c2)/2.0
     der[1,5] =  (9.0*(6.0*c1+6.0*c2-5.0)*c2)/2.0
     der[1,4] = -(27.0*c1*c1+54.0*c1*c2-36.0*c1+27.0*c2*c2-36.0*c2+11.0)/2.0
     der[1,3] =  (9.0*(9.0*c1*c1+12.0*c1*c2-10.0*c1+3.0*c2*c2-5.0*c2+2.0))/2.0
     der[1,2] = -(9.0*(9.0*c1*c1+6.0*c1*c2-8.0*c1-c2+1.0))/2.0
     der[1,10] = -27.0*(((c2-1.0)+c1)+c1)*c2
     der[2,1] = 0.0
     der[2,9] =  (9.0*(3.0*c1-1.0)*c1)/2.0
     der[2,8] =  (9.0*(6.0*c2-1.0)*c1)/2.0
     der[2,7] = (27.0*c2*c2-18.0*c2+2.0)/2.0
     der[2,6] = -(9.0*((c1+c2-1.0)*(6.0*c2-1.0)+(3.0*c2-1.0)*c2))/2.0
     der[2,5] =  (9.0*(3.0*c1*c1+12.0*c1*c2-5.0*c1+9.0*c2*c2-10.0*c2+2.0))/2.0
     der[2,4] = -(27.0*c1*c1+54.0*c1*c2-36.0*c1+27.0*c2*c2-36.0*c2+11.0)/2.0
     der[2,3] =  (9.0*(6.0*c1+6.0*c2-5.0)*c1)/2.0
     der[2,2] = -(9.0*(3.0*c1-1.0)*c1)/2.0
     der[2,10] = -27.0*(((c2-1.0)+c1)+c2)*c1
   elseif nod == 15
     t1=c1-0.25
     t2=c1-0.5
     t3=c1-0.75
     t4=c2-0.25
     t5=c2-0.5
     t6=c2-0.75
     t7=c3-0.25
     t8=c3-0.5
     t9=c3-0.75
     der[1,1] = 32.0/3.0*(t2*t3*(t1+c1)+c1*t1*(t3+t2))
     der[1,12] = 128.0/3.0*c2*(t2*(t1+c1)+c1*t1)
     der[1,11] = 64.0*c2*t4*(t1+c1)
     der[1,10] = 128.0/3.0*c2*t4*t5
     der[1,9] = 0.0
     der[1,8] = -128.0/3.0*c2*t4*t5
     der[1,7] = -64.0*c2*t4*(t7+c3)
     der[1,6] = -128.0/3.0*c2*(t8*(t7+c3)+c3*t7)
     der[1,5] = -32.0/3.0*(t8*t9*(t7+c3)+c3*t7*(t8+t9))
     der[1,4] = 128.0/3.0*(c3*t7*t8-c1*(t8*(t7+c3)+c3*t7))
     der[1,3] = 64.0*(c3*t7*(t1+c1)-c1*t1*(t7+c3))
     der[1,2] = 128.0/3.0*(c3*(t2*(t1+c1)+c1*t1)-c1*t1*t2)
     der[1,13] = 128.0*c2*(c3*(t1+c1)-c1*t1)
     der[1,15] = 128.0*c2*t4*(c3-c1)
     der[1,14] = 128.0*c2*(c3*t7-c1*(t7+c3))
     der[2,1] = 0.0
     der[2,12] = 128.0/3.0*c1*t1*t2
     der[2,11] = 64.0*c1*t1*(t4+c2)
     der[2,10] = 128.0/3.0*c1*(t5*(t4+c2)+c2*t4)
     der[2,9] = 32.0/3.0*(t5*t6*(t4+c2)+c2*t4*(t6+t5))
     der[2,8] = 128.0/3.0*((c3*(t5*(t4+c2)+c2*t4))-c2*t4*t5)
     der[2,7] = 64.0*(c3*t7*(t4+c2)-c2*t4*(t7+c3))
     der[2,6] = 128.0/3.0*(c3*t7*t8-c2*(t8*(t7+c3)+c3*t7))
     der[2,5] = -32.0/3.0*(t8*t9*(t7+c3)+c3*t7*(t8+t9))
     der[2,4] = -128.0/3.0*c1*(t8*(t7+c3)+c3*t7)
     der[2,3] = -64.0*c1*t1*(t7+c3)
     der[2,2] = -128.0/3.0*c1*t1*t2
     der[2,13] = 128.0*c1*t1*(c3-c2)
     der[2,15] = 128.0*c1*(c3*(t4+c2)-c2*t4)
     der[2,14] = 128.0*c1*(c3*t7-c2*(c3+t7))
   elseif nod == 4
     der[1,1] = -etam
     der[1,2] = -etap
     der[1,3] = etap
     der[1,4] = etam
     der[2,1] = -xim
     der[2,2] = xim
     der[2,3] = xip
     der[2,4] = -xip
   elseif nod == 5
     der[1,1] = -etam+0.5*xi*(1.0-eta*eta)
     der[1,2] = -etap+0.5*xi*(1.0-eta*eta)
     der[1,3] = etap+0.5*xi*(1.0-eta*eta)
     der[1,4] = etam+0.5*xi*(1.0-eta*eta)
     der[1,5] = -2.0*xi*(1.0-eta*eta)
     der[2,1] = -xim+0.5*eta*(1.0-xi*xi)
     der[2,2] = xim+0.5*eta*(1.0-xi*xi)
     der[2,3] = xip+0.5*eta*(1.0-xi*xi)
     der[2,4] = -xip+0.5*eta*(1.0-xi*xi)
     der[2,5] = -2.0*eta*(1.0-xi*xi)
   elseif nod == 8
     der[1,1] = etam*(2.0*xi+eta)
     der[1,2] = -8.0*etam*etap
     der[1,3] = etap*(2.0*xi-eta)
     der[1,4] = -4.0*etap*xi
     der[1,5] = etap*(2.0*xi+eta)
     der[1,6] = 8.0*etap*etam
     der[1,7] = etam*(2.0*xi-eta)
     der[1,8] = -4.0*etam*xi
     der[2,1] = xim*(xi+2.0*eta)
     der[2,2] = -4.0*xim*eta
     der[2,3] = xim*(2.0*eta-xi)
     der[2,4] = 8.0*xim*xip
     der[2,5] = xip*(xi+2.0*eta)
     der[2,6] = -4.0*xip*eta
     der[2,7] = xip*(2.0*eta-xi)
     der[2,8] = -8.0*xim*xip
   elseif nod == 9
     etam=eta-1.0
     etap=eta+1.0
     xim=xi-1.0
     xip=xi+1.0
     der[1,1] = 0.25*x2m1*eta*etam
     der[1,2] = -0.5*x2m1*etap*etam
     der[1,3] = 0.25*x2m1*eta*etap
     der[1,4] = -xi*eta*etap
     der[1,5] = 0.25*x2p1*eta*etap
     der[1,6] = -0.5*x2p1*etap*etam
     der[1,7] = 0.25*x2p1*eta*etam
     der[1,8] = -xi*eta*etam
     der[1,9] = 2.0*xi*etap*etam
     der[2,1] = 0.25*xi*xim*e2m1
     der[2,2] = -xi*xim*eta
     der[2,3] = 0.25*xi*xim*e2p1
     der[2,4] = -0.5*xip*xim*e2p1
     der[2,5] = 0.25*xi*xip*e2p1
     der[2,6] = -xi*xip*eta
     der[2,7] = 0.25*xi*xip*e2m1
     der[2,8] = -0.5*xip*xim*e2m1
     der[2,9] = 2.0*xip*xim*eta
   else
     println("Wrong number of nodes in shape_der.")
   end
  elseif ndim == 3                                          # 3.0 dimensional fin_els
   xi=points[i, 1]
   eta=points[i, 2]
   zeta=points[i, 3]
   etam=1.0-eta
   xim=1.0-xi
   zetam=1.0-zeta
   etap=eta+1.0
   xip=xi+1.0
   zetap=zeta+1.0
   if nod == 4
     der[1:3,1:4] .= 0.0
     der[1,1] = 1.0
     der[2,2] = 1.0
     der[3,3] = 1.0
     der[1,4] = -1.0
     der[2,4] = -1.0
     der[3,4] = -1.0
   elseif nod == 8
     der[1,1] = -0.125*etam*zetam
     der[1,2] = -0.125*etam*zetap
     der[1,3] =  0.125*etam*zetap
     der[1,4] =  0.125*etam*zetam
     der[1,5] = -0.125*etap*zetam
     der[1,6] = -0.125*etap*zetap
     der[1,7] =  0.125*etap*zetap
     der[1,8] =  0.125*etap*zetam
     der[2,1] = -0.125*xim*zetam
     der[2,2] = -0.125*xim*zetap
     der[2,3] = -0.125*xip*zetap
     der[2,4] = -0.125*xip*zetam
     der[2,5] =  0.125*xim*zetam
     der[2,6] =  0.125*xim*zetap
     der[2,7] =  0.125*xip*zetap
     der[2,8] =  0.125*xip*zetam
     der[3,1] = -0.125*xim*etam
     der[3,2] =  0.125*xim*etam
     der[3,3] =  0.125*xip*etam
     der[3,4] = -0.125*xip*etam
     der[3,5] = -0.125*xim*etap
     der[3,6] =  0.125*xim*etap
     der[3,7] =  0.125*xip*etap
     der[3,8] = -0.125*xip*etap
   elseif nod == 14                                     # type 6 fin_el
     der[1,1] =  (2.0*xi*eta+2.0*xi*zeta+4.0*xi+eta*zeta+eta+zeta)*
       (eta-1.0)*(zeta-1.0)/8.0
     der[1,2] = -(2.0*xi*eta-2.0*xi*zeta+4.0*xi-eta*zeta+eta-zeta)*
       (eta-1.0)*(zeta+1.0)/8.0
     der[1,3] = -(2.0*xi*eta-2.0*xi*zeta+4.0*xi+eta*zeta-eta+zeta)*
       (eta-1.0)*(zeta+1.0)/8.0
     der[1,4] =  (2.0*xi*eta+2.0*xi*zeta+4.0*xi-eta*zeta-eta-zeta)*
       (eta-1.0)*(zeta-1.0)/8.0
     der[1,5] =  -(eta-1.0)*(zeta+1.0)*(zeta-1.0)*xi
     der[1,6] = -(eta+1.0)*(eta-1.0)*(zeta+1.0)*(zeta-1.0)/2.0
     der[1,7] =   (eta+1.0)*(eta-1.0)*(zeta+1.0)*xi
     der[1,8] =  (eta+1.0)*(eta-1.0)*(zeta+1.0)*(zeta-1.0)/2.0
     der[1,9] =  -(eta+1.0)*(eta-1.0)*(zeta-1.0)*xi
     der[1,10] =  (2.0*xi*eta-2.0*xi*zeta-4.0*xi+eta*zeta+eta-zeta)*
       (eta+1.0)*(zeta-1.0)/8.0
     der[1,11] = -(2.0*xi*eta+2.0*xi*zeta-4.0*xi-eta*zeta+eta+zeta)*
       (eta+1.0)*(zeta+1.0)/8.0
     der[1,12] = -(2.0*xi*eta+2.0*xi*zeta-4.0*xi+eta*zeta-eta-zeta)*
       (eta+1.0)*(zeta+1.0)/8.0
     der[1,13] =  (2.0*xi*eta-2.0*xi*zeta-4.0*xi-eta*zeta-eta+zeta)*
       (eta+1.0)*(zeta-1.0)/8.0
     der[1,14] =   (eta+1.0)*(zeta+1.0)*(zeta-1.0)*xi
     der[2,1] =  (2.0*xi*eta+xi*zeta+xi+2.0*eta*zeta+4.0*eta+zeta)*
       (xi-1.0)*(zeta-1.0)/8.0
     der[2,2] = -(2.0*xi*eta-xi*zeta+xi-2.0*eta*zeta+4.0*eta-zeta)*
       (xi-1.0)*(zeta+1.0)/8.0
     der[2,3] = -(2.0*xi*eta-xi*zeta+xi+2.0*eta*zeta-4.0*eta+zeta)*
       (xi+1.0)*(zeta+1.0)/8.0
     der[2,4] =  (2.0*xi*eta+xi*zeta+xi-2.0*eta*zeta-4.0*eta-zeta)*
       (xi+1.0)*(zeta-1.0)/8.0
     der[2,5] = -(xi+1.0)*(xi-1.0)*(zeta+1.0)*(zeta-1.0)/2.0
     der[2,6] =  -(xi-1.0)*(zeta+1.0)*(zeta-1.0)*eta
     der[2,7] =   (xi+1.0)*(xi-1.0)*(zeta+1.0)*eta
     der[2,8] =   (xi+1.0)*(zeta+1.0)*(zeta-1.0)*eta
     der[2,9] =  -(xi+1.0)*(xi-1.0)*(zeta-1.0)*eta
     der[2,10] =  (2.0*xi*eta-xi*zeta-xi+2.0*eta*zeta+4.0*eta-zeta)*
       (xi-1.0)*(zeta-1.0)/8.0
     der[2,11] = -(2.0*xi*eta+xi*zeta-xi-2.0*eta*zeta+4.0*eta+zeta)*
       (xi-1.0)*(zeta+1.0)/8.0
     der[2,12] = -(2.0*xi*eta+xi*zeta-xi+2.0*eta*zeta-4.0*eta-zeta)*
       (xi+1.0)*(zeta+1.0)/8.0
     der[2,13] =  (2.0*xi*eta-xi*zeta-xi-2.0*eta*zeta-4.0*eta+zeta)*
       (xi+1.0)*(zeta-1.0)/8.0
     der[2,14] =  (xi+1.0)*(xi-1.0)*(zeta+1.0)*(zeta-1.0)/2.0
     der[3,1] =  (xi*eta+2.0*xi*zeta+xi+2.0*eta*zeta+eta+4.0*zeta)*
       (xi-1.0)*(eta-1.0)/8.0
     der[3,2] = -(xi*eta-2.0*xi*zeta+xi-2.0*eta*zeta+eta-4.0*zeta)*
       (xi-1.0)*(eta-1.0)/8.0
     der[3,3] = -(xi*eta-2.0*xi*zeta+xi+2.0*eta*zeta-eta+4.0*zeta)*
       (xi+1.0)*(eta-1.0)/8.0
     der[3,4] =  (xi*eta+2.0*xi*zeta+xi-2.0*eta*zeta-eta-4.0*zeta)*
       (xi+1.0)*(eta-1.0)/8.0
     der[3,5] =  -(xi+1.0)*(xi-1.0)*(eta-1.0)*zeta
     der[3,6] =  -(xi-1.0)*(eta+1.0)*(eta-1.0)*zeta
     der[3,7] =  (xi+1.0)*(xi-1.0)*(eta+1.0)*(eta-1.0)/2.0
     der[3,8] =   (xi+1.0)*(eta+1.0)*(eta-1.0)*zeta
     der[3,9] = -(xi+1.0)*(xi-1.0)*(eta+1.0)*(eta-1.0)/2.0
     der[3,10] =  (xi*eta-2.0*xi*zeta-xi+2.0*eta*zeta+eta-4.0*zeta)*
       (xi-1.0)*(eta+1.0)/8.0
     der[3,11] = -(xi*eta+2.0*xi*zeta-xi-2.0*eta*zeta+eta+4.0*zeta)*
       (xi-1.0)*(eta+1.0)/8.0
     der[3,12] = -(xi*eta+2.0*xi*zeta-xi+2.0*eta*zeta-eta-4.0*zeta)*
       (xi+1.0)*(eta+1.0)/8.0
     der[3,13] =  (xi*eta-2.0*xi*zeta-xi-2.0*eta*zeta-eta+4.0*zeta)*
       (xi+1.0)*(eta+1.0)/8.0
     der[3,14] =   (xi+1.0)*(xi-1.0)*(eta+1.0)*zeta
   elseif nod == 20
     xii = [-1,-1,-1,0,1,1,1,0,-1,-1,1,1,-1,-1,-1,0,1,1,1,0]
     etai = [-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,1,1,1,1,1,1,1,1]
     zetai = [-1,0,1,1,1,0,-1,-1,-1,1,1,-1,-1,0,1,1,1,0,-1,-1]
     for l in 1:20
       xi0=xi*xii[l]
       eta0=eta*etai[l]
       zeta0=zeta*zetai[l]
       if l == 4 || l == 8 || l == 16 || l == 20
         der[1,l] = -0.5*xi*(1.0+eta0)*(1.0+zeta0)
         der[2,l] = 0.25*etai[l]*(1.0-xi*xi)*(1.0+zeta0)
         der[3,l] = 0.25*zetai[l]*(1.0-xi*xi)*(1.0+eta0)
       elseif l >= 9 && l <= 12
         der[1,l] = 0.25*xii[l]*(1.0-eta*eta)*(1.0+zeta0)
         der[2,l] = -0.5*eta*(1.0+xi0)*(1.0+zeta0)
         der[3,l] = 0.25*zetai[l]*(1.0+xi0)*(1.0-eta*eta)
       elseif l == 2 || l == 6 || l == 14 || l == 18
         der[1,l] = 0.25*xii[l]*(1.0+eta0)*(1.0-zeta*zeta)
         der[2,l] = 0.25*etai[l]*(1.0+xi0)*(1.0-zeta*zeta)
         der[3,l] = -0.5*zeta*(1.0+xi0)*(1.0+eta0)
       else
         der[1,l] = 0.125*xii[l]*(1.0+eta0)*(1.0+zeta0)*
           (2.0*xi0+eta0+zeta0-1.0)
         der[2,l] = 0.125*etai[l]*(1.0+xi0)*(1.0+zeta0)*
           (xi0+2.0*eta0+zeta0-1.0)
         der[3,l] = 0.125*zetai[l]*(1.0+xi0)*(1.0+eta0)*
           (xi0+eta0+2.0*zeta0-1.0)
       end
     end
   else
     println("Wrong number of nodes in shape_der.")
   end
  else
   println("Wrong number of dimensions in shape_der.")
  end
end


text:
using SeisMain
using Test

println("Testing test_bin.jl")

download("http://seismic.physics.ualberta.ca/data/prestack_section.su","section.su");

SegyToSeis("section.su","section",format="su",input_type="ieee");

param1 = Dict( :dmx=>15, :dmy=>15, :dh=>30, :daz=>45 );

param2 = Dict(:style=>"mxmyhaz", :min_imx=>10,:max_imx=>100, :min_imy=>35, :max_imy=>45,
               :min_ih=>1, :max_ih=>6, :min_iaz=>0, :max_iaz=>7);

SeisGeometry("section"; param1...)

SeisBinHeaders("section","section_bin"; param1..., param2...);

SeisBinData("section","section_bin"; param1..., param2...);

db,hb,eb=SeisRead("section_bin");

N =size(db)

@test N == (251, 91, 11, 6, 8)


text:
###############################################################################
#
#   Float.jl : Additional Nemo functionality for Julia AbstractFloat types
#
###############################################################################

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

JuliaRealField = Floats{BigFloat}()

RDF = Floats{Float64}()

parent(a::T) where T <: AbstractFloat = Floats{T}()

elem_type(::Type{Floats{T}}) where T <: AbstractFloat = T

parent_type(::Type{T}) where T <: AbstractFloat = Floats{T}

base_ring(a::AbstractFloat) = Union{}

isdomain_type(::Type{T}) where T <: AbstractFloat = true

isexact_type(::Type{T}) where T <: AbstractFloat = false

###############################################################################
#
#   Basic manipulation
#
###############################################################################

zero(::Floats{T}) where T <: AbstractFloat = T(0)

one(::Floats{T}) where T <: AbstractFloat = T(1)

if VERSION < v"0.7.0-DEV.1144"
isone(a::AbstractFloat) = a == 1
end

isunit(a::AbstractFloat) = a != 0

canonical_unit(a::AbstractFloat) = a

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, R::Floats)
   print(io, "Floats")
end

needs_parentheses(::AbstractFloat) = false

isnegative(a::AbstractFloat) = a < 0

show_minus_one(::Type{T}) where T <: AbstractFloat = false

###############################################################################
#
#   Exact division
#
###############################################################################

divexact(a::T, b::T) where T <: AbstractFloat = a/b

divexact(a::AbstractFloat, b::Union{Signed, Unsigned}) = a/b

divexact(a::BigFloat, b::BigInt) = a/b

divexact(a::Union{Signed, Unsigned}, b::AbstractFloat) = a/b

divexact(a::BigInt, b::BigFloat) = a/b

divexact(a::AbstractFloat, b::Rational{T}) where T <: Union{Signed, Unsigned} = a/b

divexact(a::BigFloat, b::Rational{BigInt}) = a/b

divexact(a::Rational{T}, b::AbstractFloat) where T <: Union{Signed, Unsigned} = a/b

divexact(a::Rational{BigInt}, b::BigFloat) = a/b

###############################################################################
#
#   GCD
#
###############################################################################

function gcd(a::T, b::T) where T <: AbstractFloat
   if a == 0 && b == 0
      return T(0)
   else
      return T(1)
   end
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(a::AbstractFloat)
   return 0
end

function zero!(a::BigFloat)
   ccall((:mpfr_set_si, :libmpfr), Void,
         (Ref{BigFloat}, Int, Int32), a, 0, Base.MPFR.ROUNDING_MODE[])
   return a
end

function mul!(a::T, b::T, c::T) where T <: AbstractFloat
   return b*c
end

function mul!(a::BigFloat, b::BigFloat, c::BigFloat)
   ccall((:mpfr_mul, :libmpfr), Void,
         (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32),
                 a, b, c, Base.MPFR.ROUNDING_MODE[])
   return a
end

function add!(a::T, b::T, c::T) where T <: AbstractFloat
   return b + c
end

function add!(a::BigFloat, b::BigFloat, c::BigFloat)
   ccall((:mpfr_add, :libmpfr), Void,
         (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32),
                 a, b, c, Base.MPFR.ROUNDING_MODE[])
   return a
end

function addeq!(a::T, b::T) where T <: AbstractFloat
   return a + b
end

function addeq!(a::BigFloat, b::BigFloat)
   ccall((:mpfr_add, :libmpfr), Void,
         (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32),
                 a, a, b, Base.MPFR.ROUNDING_MODE[])
   return a
end

function addmul!(a::T, b::T, c::T, d::T) where T <: AbstractFloat
   return a + b*c
end

function addmul!(a::BigFloat, b::BigFloat, c::BigFloat, d::BigFloat)
   ccall((:mpfr_fma, :libmpfr), Void,
         (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32),
                 a, b, c, a, Base.MPFR.ROUNDING_MODE[])
   return a
end

function addmul!(a::T, b::T, c::T) where T <: AbstractFloat # special case, no temporary required
   return a + b*c
end

function addmul!(a::BigFloat, b::BigFloat, c::BigFloat) # special case, no temporary required
   ccall((:mpfr_fma, :libmpfr), Void,
         (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32),
                 a, b, c, a, Base.MPFR.ROUNDING_MODE[])
   return a
end

###############################################################################
#
#   Random generation
#
###############################################################################

function rand(R::Floats, n::UnitRange{AbstractFloat})
   return R(n.start + rand(Float64)*(n.stop - n.start))
end

function rand(R::Floats, n::UnitRange{Int})
   return R(n.start + rand(Float64)*(n.stop - n.start))
end

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::Floats{T})() where T <: AbstractFloat
   return T(0)
end

function (a::Floats{T})(b::AbstractFloat) where T <: AbstractFloat
   return T(b)
end

function (a::Floats{T})(b::Union{Signed, Unsigned}) where T <: AbstractFloat
   return T(b)
end

function (a::Floats{BigFloat})(b::BigInt)
   return BigFloat(b)
end

function (a::Floats{T})(b::Rational{S}) where {S <: Union{Signed, Unsigned}, T <: AbstractFloat}
   return T(b)
end

function (a::Floats{BigFloat})(b::Rational{BigInt})
   return BigFloat(b)
end


text:
using Test
using TightBindingLattice

@testset "Permutation" begin
    @test_throws ArgumentError Permutation([1,2,4])
    @test_throws OverflowError Permutation([mod(x, 4096)+1 for x in 1:4096])
    p0 = Permutation([1,2,3,4])
    p1 = Permutation([2,3,4,1])
    p2 = Permutation([3,4,1,2])
    p3 = Permutation([4,1,2,3])

    @test p1 * p2 == p3
    @test p1 != p3
    @test p1^0 == p0
    @test p1^1 == p1
    @test p1^2 == p2
    @test p1^3 == p3

    @test p0.order == 1
    @test p1.order == 4
    @test p2.order == 2
    @test p3.order == 4

    @test inv(p0) == p0
    @test inv(p1) == p3
    @test inv(p2) == p2
    @test inv(p3) == p1

    @test p0 == Permutation([1,2,3,4])
    @test isless(p0, p1)
    @test isless(p1, p3)

    @test p0(1) == 1 && p0(2) == 2 && p0(3) == 3 && p0(4) == 4
    @test p1(1) == 2 && p1(2) == 3 && p1(3) == 4 && p1(4) == 1

    @test_throws ArgumentError Permutation([1,2,3,4]) * Permutation([1,2,3,4,5])
    @test hash(Permutation(Int[1,2,3,4])) != hash(Int[1,2,3,4])

    @test generate_group(p1) == Set([p0, p1, p2, p3])
    @test generate_group(p2) == Set([p0, p2])
    @test generate_group(p1, p2) == Set([p0, p1, p2, p3])
end


text:
using Documenter
using Brownian

makedocs(
    sitename = "Brownian",
    format = Documenter.HTML(),
    modules = [Brownian]
)

deploydocs(
    repo = "github.com/pnavaro/Brownian.jl.git"
)


text:
@testset "Cross mapping" begin

	@testset "Prediction lags" begin
	    x, y = rand(100), rand(100)
	    crossmap(x, y, η = 1)
	    crossmap(x, y, η = -1)
	    crossmap(x, y, η = 5)
	    crossmap(x, y, η = -5)
	end

	@testset "Embedding params" begin
	    x, y = rand(100), rand(100)
	    @test_throws DomainError crossmap(x, y, dim = -3)
	    @test_throws DomainError crossmap(x, y, dim = 100, τ = 2)
	    crossmap(x, y, dim = 3)
	    crossmap(x, y, dim = 10, τ = 2)
	end

	@testset "Exclusion radii" begin
	    x, y = rand(100), rand(100)
	    [crossmap(x, y, theiler_window = i) for i in rand(1:25, 10)]
	    #@test_throws DomainError crossmap(x, y, theiler_window = -1)
	end

	@testset "Replacements" begin
	    x, y = rand(100), rand(100)
	    crossmap(x, y, replace = false)
	    crossmap(x, y, replace = true)
	end

	@testset "Correspondence measures" begin
	    @testset "$i" for i in 1:100
	        x, y = rand(100), rand(100)
	        @test all(crossmap(x, y, correspondence_measure = StatsBase.rmsd) .>= 0)
	        @test all([-1 <= x <= 1 for x in crossmap(x, y, correspondence_measure = StatsBase.cor)])
	    end
	end

end


text:
import Ghost: Tape, V, inputs!, rebind!, mkcall, primitivize!

@testset "tape" begin
    # rebind!
    tape = Tape()
    _, v1, v2 = inputs!(tape, nothing, 3.0, 5.0)
    v3 = push!(tape, mkcall(*, v1, 2))
    st = Dict(v1.id => v2.id)
    rebind!(tape, st)
    @test tape[v3].args[1].id == v2.id

    # variable equality
    @test tape[v3].args[1] == v2         # bound var
    @test tape[v3].args[1] != V(v2.id)   # unbound var

    # mkcall value calculation
    c = mkcall(*, 2.0, v1)                # bound var
    @test c.val == 2.0 * tape[v1].val
    c = mkcall(*, 2.0, V(100))            # unbound var
    @test c.val === missing
    c = mkcall(*, 2.0, V(100); val=10.0)  # manual value
    @test c.val == 10.0

    # push!, insert!; var hash
    tape = Tape()
    a1, a2, a3 = inputs!(tape, nothing, 2.0, 5.0)
    r = push!(tape, mkcall(*, a2, a3))
    @test tape[r].val == 10.0

    dct = Dict(r => :r)
    ops = [mkcall(+, a2, 1), mkcall(+, a3, 1)]
    v1, v2 = insert!(tape, 4, ops...)
    @test r.id == 6
    @test dct[r] == :r

    tape[r] = mkcall(*, v1, v2)
    @test tape[r].val == 18.0

    v2.id = 100
    @test tape[r].args[2].id == 100

    # replace!
    tape2 = deepcopy(tape)
    tape3 = deepcopy(tape)

    op1 = mkcall(*, V(2), 2)
    op2 = mkcall(+, V(op1), 1)
    z = replace!(tape, 4 => [op1, op2]; rebind_to=2)
    @test tape[V(7)].args[1].id == op2.id
    @test z isa V

    replace!(tape2, V(4) => [op1, op2]; rebind_to=2)
    @test tape2[V(7)].args[1].id == op2.id

    replace!(tape3, tape3[V(4)] => [op1, op2]; rebind_to=2)
    @test tape3[V(7)].args[1].id == op2.id

    # deleteat!
    tape = Tape()
    _, v1, v2 = inputs!(tape, nothing, 3.0, 5.0)
    v3 = push!(tape, mkcall(*, v1, 2))
    v4 = push!(tape, mkcall(+, v3, v1))
    v5 = push!(tape, mkcall(+, v4, v2))
    v6 = push!(tape, mkcall(+, v4, v1))
    tape2, tape3 = deepcopy(tape), deepcopy(tape)

    deleteat!(tape, 5; rebind_to = 1)
    @test tape[V(5)].args[1].id == 1
    @test tape[V(6)].args[1].id == 1
    @test all(x -> x[1].id == x[2], zip(tape, 1:length(tape)))

    deleteat!(tape2, 5; rebind_to = 1)
    @test tape2[V(5)].args[1].id == 1
    @test tape2[V(6)].args[1].id == 1
    @test all(x -> x[1].id == x[2], zip(tape2, 1:length(tape2)))

    deleteat!(tape3, 5; rebind_to = 1)
    @test tape3[V(5)].args[1].id == 1
    @test tape3[V(6)].args[1].id == 1
    @test all(x -> x[1].id == x[2], zip(tape3, 1:length(tape3)))

    tape = Tape()
    _, v1, v2 = inputs!(tape, nothing, 3.0, 5.0)
    v3 = push!(tape, mkcall(println, "Test"))
    deleteat!(tape, v3)
    @test length(tape) == 3

    # primitivize!
    f(x) = 2x - 1
    g(x) = f(x) + 5

    tape = Tape()
    _, x = inputs!(tape, g, 3.0)
    y = push!(tape, mkcall(f, x))
    z = push!(tape, mkcall(+, y, 5))
    tape.result = z

    primitivize!(tape)

    @test length(tape) == 5
    @test tape[V(3)].fn == *
    @test tape[V(4)].fn == -

end


text:
abstract type Runner end

"""
Holds all the plumbing required for a running process or function
"""
struct Plumbing
    inp::Pipe                               # stdin
    out::Pipe                               # stdout
    err::Pipe                               # stderr

    function Plumbing()
        new(Pipe(), Pipe(), Pipe())
    end
end

function close(plumbing::Plumbing; inputendsonly::Bool=false)
    if inputendsonly
        close(plumbing.out.in)
        close(plumbing.err.in)
    else
        close(plumbing.inp)
        close(plumbing.out)
        close(plumbing.err)
    end
    nothing
end

function link(plumbing::Plumbing)
    Base.link_pipe!(plumbing.inp)
    Base.link_pipe!(plumbing.out)
    Base.link_pipe!(plumbing.err)
    nothing
end

"""
Runs fast CGI commands as processes. Holds the runner process and allows waiting and terminating the process.
"""
mutable struct ProcessRunner <: Runner
    process::Union{Base.Process,Nothing}    # launched process
    timeout::Int                            # time the process out after these many seconds
    timer::Union{Timer,Nothing}             # timer to kill the running process after timeout

    function ProcessRunner(timeout::Int=-1)
        new(nothing, timeout, nothing)
    end
end

function killproc(runner::ProcessRunner)
    if runner.process !== nothing
        try
            kill(runner.process)
        catch ex
            @warn("failed to kill process with SIGTERM, will issue SIGKILL")
            try
                kill(runner.process, Base.SIGKILL)
            catch ex
                @warn("failed to kill process with SIGKILL, abandoning")
            end
        end
        runner.process = nothing
    end
    closetimer(runner)
    nothing
end

function waitproc(runner::ProcessRunner)
    exitcode = 408
    try
        wait(runner.process)
        if runner.process !== nothing
            exitcode = runner.process.exitcode
        end
    catch ex
        @warn("exception waiting for process", ex)
    end
    runner.process = nothing
    closetimer(runner)
    exitcode
end

function launchproc(runner::ProcessRunner, plumbing::Plumbing, params::Dict{String,String})
    cmdpath = getcommand(params)
    @debug("launching command", cmdpath)
    if !isfile(cmdpath)
        err = "cmdpath not found: $cmdpath"
        @warn(err)
        return 404, err
    end
    try
        cmd = Cmd(`$cmdpath`; env=params)
        runner.process = run(pipeline(cmd, stdin=plumbing.inp, stdout=plumbing.out, stderr=plumbing.err), wait=false)
        setuptimer(runner, plumbing.out, plumbing.err)
        return 0, ""
    catch ex
        @warn("process exception", cmdpath, params, ex)
        return 500, "process exception"
    end
end

"""
Runs fast CGI commands as Julia functions. Holds the runner task and allows waiting and terminating (interrupting) the task.
"""
mutable struct FunctionRunner <: Runner
    process::Union{Task,Nothing}     # launched process
    timeout::Int                     # time the process out after these many seconds
    timer::Union{Timer,Nothing}      # timer to kill the running process after timeout

    function FunctionRunner(timeout::Int=-1)
        new(nothing, timeout, nothing)
    end
end

function killproc(runner::FunctionRunner)
    if (runner.process !== nothing) && !istaskdone(runner.process)
        try
            Base.throwto(runner.process, InterruptException())
        catch ex
            # ignore
        end
    end
    runner.process = nothing
    closetimer(runner)
    nothing
end

function waitproc(runner::FunctionRunner)
    result = 408
    try
        wait(runner.process)
        if runner.process !== nothing
            result = Int(fetch(runner.process))
        end
    catch ex
        # ignore (may have been interrupted)
        @warn("exception waiting for process", ex)
    end
    runner.process = nothing
    closetimer(runner)
    result
end

function launchproc(runner::FunctionRunner, plumbing::Plumbing, params::Dict{String,String})
    cmdpath = getcommand(params)
    @debug("launching command", cmdpath)

    cmd = try
        T = Main
        for t in split(cmdpath, ".")
            T = Base.eval(T, Symbol(t))
        end
        T::Function
    catch ex
        err = "Can not resolve $cmdpath"
        @warn(err, ex)
        return 404, err
    end
    try
        link(plumbing)
        runner.process = @async cmd(params, Base.pipe_reader(plumbing.inp), Base.pipe_writer(plumbing.out), Base.pipe_writer(plumbing.err))
        setuptimer(runner, Base.pipe_writer(plumbing.out), Base.pipe_writer(plumbing.err))
        return 0, ""
    catch ex
        @warn("process exception", cmdpath, params, ex)
        return 500, "process exception"
    end
end

function closetimer(runner::Union{ProcessRunner,FunctionRunner})
    if runner.timer !== nothing
        try
            close(runner.timer)
        catch ex
            @warn("exception closing timer", ex)
        end
        runner.timer = nothing
    end
    nothing
end

function setuptimer(runner::Union{ProcessRunner,FunctionRunner}, out, err)
    if runner.timeout > 0
        runner.timer = Timer(runner.timeout) do timer
            close(out)
            close(err)
            killproc(runner)
        end
    end
    nothing
end


text:
function calc_fullevo(initialstate,time_evo_array)
  state_prob = zeros(Complex128,length(time_evo_array), size(initialstate, 1),size(initialstate, 2))
  # Time evolve the system density matrix and "measure" the qubit at each time sample taken. This is done by
  # tracing out the resonator and keeping the excited state entry of the density matrix (ie, entry 2,2 in the matrix)
  for i=1:length(time_evo_array)
    state_prob[i, :, :] = partialtrace(time_evo_array[i] * initialstate * time_evo_array[i]',size(initialstate, 1))
  end
  return state_prob
end


text:
@doc raw"""
    blacken(white_rbm::WhiteRBM)

Constructs a plain `RBM` equivalent to the given `white_rbm`.
The energies assigned by the two models differ by a constant amount,

```math
E(v,h) - \tilde{E}(v,h) = \mathbf{a}^\top \mathbb{A}^\top \tilde{\mathbb{W}} \mathbb{B}\mathbf{b}
```

where ``\tilde{E}(v,h)`` is the energy assigned by `white_rbm` and ``E(v,h)`` is the energy
assigned by the `RBM` constructed by this method.

This is the inverse operation of [`whiten`](@ref), always returning an ordinary `RBM`.

To construct an `RBM` that simply neglects the transformations,
call `RBM(white_rbm)` instead.
"""
blacken(white_rbm::WhiteRBM) = RBM(whiten(white_rbm))
blacken(rbm::RBM) = rbm

@doc raw"""
    whiten(rbm::RBM, affine_v, affine_h)

Constructs a `WhiteRBM` equivalent to the given `rbm`.
The energies assigned by the two models differ by a constant amount,

```math
E(v,h) - \tilde{E}(v,h) = \mathbf{a}^\top \mathbb{A}^\top \tilde{\mathbb{W}} \mathbb{B}\mathbf{b}
```

where ``E(v,h)`` is the energy assigned by the original `rbm`, and
``\tilde{E}(v,h)`` is the energy assigned by the returned `WhiteRBM`.

This is the inverse operation of [`blacken`](@ref), always returning a `WhiteRBM`.

To construct a `WhiteRBM` that simply includes these transforms,
call `WhiteRBM(rbm, affine_v, affine_h)` instead.
"""
whiten(rbm::RBM, affine_v::Affine, affine_h::Affine) = whiten(whiten(rbm), affine_v, affine_h)
whiten(rbm::RBM) = WhiteRBM(rbm)

function whiten(white_rbm::AffineRBM, affine_v::Affine, affine_h::Affine)
    @assert length(visible(white_rbm)) == length(affine_v.u)
    @assert length(hidden(white_rbm)) == length(affine_h.u)

    w1 = reshape(weights(white_rbm), length(visible(white_rbm)), length(hidden(white_rbm)))
    w2 = affine_v.A' \ white_rbm.affine_v.A' * w1 * white_rbm.affine_h.A / affine_h.A
    Δg = inputs_h_to_v(white_rbm, reshape(affine_h.u, size(hidden(white_rbm))))
    Δθ = inputs_v_to_h(white_rbm, reshape(affine_v.u, size(visible(white_rbm))))

    vis = shift_fields(visible(white_rbm), reshape(Δg, size(visible(white_rbm))))
    hid = shift_fields(hidden(white_rbm), reshape(Δθ, size(hidden(white_rbm))))
    rbm = RBM(vis, hid, reshape(w2, size(weights(white_rbm))))

    return WhiteRBM(rbm, affine_v, affine_h)
end

function whiten(white_rbm::WhiteRBM)
    affine_v = one(white_rbm.affine_v)
    affine_h = one(white_rbm.affine_h)
    return whiten(white_rbm, affine_v, affine_h)
end

function whiten_visible(white_rbm::WhiteRBM, affine_v::Affine)
    @assert length(visible(white_rbm)) == length(affine_v.u)

    w1 = reshape(weights(white_rbm), length(visible(white_rbm)), length(hidden(white_rbm)))
    w2 = affine_v.A' \ white_rbm.affine_v.A' * w1
    Δθ = inputs_v_to_h(white_rbm, reshape(affine_v.u, size(visible(white_rbm))))

    hid = shift_fields(hidden(white_rbm), Δθ)
    rbm = RBM(visible(white_rbm), hid, reshape(w2, size(weights(white_rbm))))

    return WhiteRBM(rbm, affine_v, white_rbm.affine_h)
end

function whiten_hidden(white_rbm::WhiteRBM, affine_h::Affine)
    @assert length(hidden(white_rbm)) == length(affine_h.u)

    w1 = reshape(weights(white_rbm), length(visible(white_rbm)), length(hidden(white_rbm)))
    w2 = w1 * white_rbm.affine_h.A / affine_h.A
    Δg = inputs_h_to_v(white_rbm, reshape(affine_h.u, size(hidden(white_rbm))))

    vis = shift_fields(visible(white_rbm), Δg)
    rbm = RBM(vis, hidden(white_rbm), reshape(w2, size(weights(white_rbm))))

    return WhiteRBM(rbm, white_rbm.affine_v, affine_h)
end

whiten_visible(rbm::RBM, affine_v::Affine) = whiten_visible(whiten(rbm), affine_v)
whiten_hidden(rbm::RBM, affine_h::Affine) = whiten_hidden(whiten(rbm), affine_h)

# In-place versions

function whiten!(rbm::AffineRBM, affine_v::Affine, affine_h::Affine)
    whiten_visible!(rbm, affine_v)
    whiten_hidden!(rbm, affine_h)
end

whiten!(rbm::WhiteRBM) = whiten!(rbm, one(rbm.affine_v), one(rbm.affine_h))

function whiten_visible!(rbm::WhiteRBM, affine_v::Affine)
    @assert length(visible(rbm)) == length(affine_v.u)
    Δθ = inputs_v_to_h(rbm, reshape(affine_v.u, size(visible(rbm))))
    shift_fields!(hidden(rbm), Δθ)
    w = reshape(weights(rbm), length(visible(rbm)), length(hidden(rbm)))
    copyto!(weights(rbm), affine_v.A' \ rbm.affine_v.A' * w)
    copy!(rbm.affine_v, affine_v)
    return rbm
end

function whiten_hidden!(rbm::WhiteRBM, affine_h::Affine)
    @assert length(hidden(rbm)) == length(affine_h.u)
    Δg = inputs_h_to_v(rbm, reshape(affine_h.u, size(hidden(rbm))))
    shift_fields!(visible(rbm), Δg)
    w = reshape(weights(rbm), length(visible(rbm)), length(hidden(rbm)))
    copyto!(weights(rbm), w * rbm.affine_h.A / affine_h.A)
    copy!(rbm.affine_h, affine_h)
    return rbm
end

"""
    energy_shift(rbm, affine_v, affine_h)

Computes the constant energy shift if the affine transformations were updated as given.
"""
function energy_shift(rbm::WhiteRBM, affine_v::Affine, affine_h::Affine)
    @assert length(visible(rbm)) == length(affine_v.u)
    @assert length(hidden(rbm)) == length(affine_h.u)
    E1 = RBMs.interaction_energy(RBM(rbm), rbm.affine_v.A * rbm.affine_v.u, rbm.affine_h.A * rbm.affine_h.u)
    E2 = RBMs.interaction_energy(RBM(rbm), rbm.affine_v.A * affine_v.u, rbm.affine_h.A * affine_h.u)
    return E2 - E1
end

function energy_shift_visible(rbm::WhiteRBM, affine_v::Affine)
    Δa = affine_v.u - rbm.affine_v.u
    return RBMs.interaction_energy(RBM(rbm), rbm.affine_v.A * Δa, rbm.affine_h.A * rbm.affine_h.u)
end

function energy_shift_hidden(rbm::WhiteRBM, affine_h::Affine)
    Δb = affine_h.u - rbm.affine_h.u
    return RBMs.interaction_energy(RBM(rbm), rbm.affine_v.A * rbm.affine_v.u, rbm.affine_h.A * Δb)
end

energy_shift(rbm::RBM, affine_v::Affine, affine_h::Affine) = energy_shift(whiten(rbm), affine_v, affine_h)
energy_shift_visible(rbm::RBM, affine_v::Affine) = energy_shift_visible(whiten(rbm), affine_v)
energy_shift_hidden(rbm::RBM, affine_h::Affine) = energy_shift_hidden(whiten(rbm), affine_h)

# Safe versions (preserve argument type)

"""
    safe_whiten(white_rbm, affine_v, affine_h)

Like `whiten(white_rbm, affine_v, affine_h)`, but ensures that the returned `WhiteRBM`
is of the same type as `white_rbm` (in particular, preserving the types of the
affine transformations).
"""
function safe_whiten(white_rbm::AffineRBM{Av,Ah}, affine_v::Av, affine_h::Ah) where {Av,Ah}
    return oftype(white_rbm, whiten(white_rbm, affine_v, affine_h))
end

safe_whiten(white_rbm::WhiteRBM) = oftype(white_rbm, whiten(white_rbm))

function safe_whiten_visible(white_rbm::AffineRBM{Av,Ah}, affine_v::Av) where {Av,Ah}
    return oftype(white_rbm, whiten_visible(white_rbm, affine_v))
end

function safe_whiten_hidden(white_rbm::AffineRBM{Av,Ah}, affine_h::Ah) where {Av,Ah}
    return oftype(white_rbm, whiten_hidden(white_rbm, affine_h))
end

# for RBM arguments fallback to normal methods
safe_whiten(rbm::RBM, affine_v::Affine, affine_h::Affine) = whiten(rbm, affine_v, affine_h)
safe_whiten(rbm::RBM) = whiten(rbm)
safe_whiten_visible(rbm::RBM, affine::Affine) = whiten_visible(rbm, affine)
safe_whiten_hidden(rbm::RBM, affine::Affine) = whiten_hidden(rbm, affine)


text:
using BinDeps
using Compat

@BinDeps.setup

const libOSXUVC_version = "master"

ignore_paths = split(strip(get(ENV, "LIBOSXUVCJL_LIBRARY_IGNORE_PATH", "")), ':')

validate = function(libpath, handle)
    for path in ignore_paths
        isempty(path) && continue
        ismatch(Regex("^$(path)"), libpath) && return false
    end
    return true
end

libOSXUVC = library_dependency("libOSXUVC", validate=validate)

# Yes, this is closed for now
github_root = "http://github.team-lab.local/ryuyamamoto/libOSXUVC"

provides(Sources,
         URI("$(github_root)/archive/$(libOSXUVC_version).tar.gz"),
         libOSXUVC,
         unpacked_dir="libOSXUVC-$(libOSXUVC_version)")

prefix = joinpath(BinDeps.depsdir(libOSXUVC), "usr")
srcdir = joinpath(BinDeps.depsdir(libOSXUVC), "src", "libOSXUVC-$(libOSXUVC_version)")

cmake_options = [
    "-DCMAKE_INSTALL_PREFIX=$prefix",
    "-DCMAKE_BUILD_TYPE=RELEASE",
]

provides(SimpleBuild,
          (@build_steps begin
              GetSources(libOSXUVC)
              @build_steps begin
                  ChangeDirectory(srcdir)
                  `mkdir -p build`
                  @build_steps begin
                      ChangeDirectory(joinpath(srcdir, "build"))
                      `rm -f CMakeCache.txt`
                      `cmake $cmake_options ..`
                      `make -j4`
                      `make install`
                  end
                end
          end), libOSXUVC, os = :Unix)

@BinDeps.install @compat Dict(:libOSXUVC => :libOSXUVC)


text:
include("imports.jl")
include("hausdorff.jl")
include("dice.jl")


text:
function makeplots(picdir)
    PyPlot.clf()
    plotting_multiscene(Plotter=PyPlot)
    fname=joinpath(picdir,"plotting_multiscene.svg")
    PyPlot.savefig(fname)
    @assert isfile(fname)
    @show "multiscene"

    PyPlot.clf()
    plotting_func3d(Plotter=PyPlot)
    fname=joinpath(picdir,"plotting_func3d.svg")
    PyPlot.savefig(fname)
    @assert isfile(fname)
    @show "func3d"

    PyPlot.clf()
    plotting_grid3d(Plotter=PyPlot)
    fname=joinpath(picdir,"plotting_grid3d.svg")
    PyPlot.savefig(fname)
    @assert isfile(fname)
    @show "grid3d"
    true
end


text:
using Clang.Generators
using ViennaRNA_jll

cd(@__DIR__)

include_dir = normpath(ViennaRNA_jll.artifact_dir, "include")
viennarna_dir = joinpath(include_dir, "ViennaRNA")

options = load_options(joinpath(@__DIR__, "generator.toml"))

args = get_default_args()
append!(args, [
    "-I$include_dir",
    "-DVRNA_DISABLE_C11_FEATURES",
])

accept_headers = [
    "equilibrium_probs.h",
    "eval.h",
    "fold_compound.h",
    "heat_capacity.h",
    "inverse.h",
    "landscape/neighbor.h",
    "MEA.h",
    "mfe.h",
    "params/io.h",
    "part_func.h",
    "plotting/layouts.h",
    "subopt.h",
    "subopt_zuker.h",
    "treedist.h",
    "utils/basic.h",
]
headers = [joinpath(viennarna_dir, header) for header in accept_headers]
# there is also an experimental `detect_headers` function for auto-detecting top-level headers in the directory
# headers = detect_headers(clang_dir, args)

# create context
ctx = create_context(headers, args, options)

# run generator
build!(ctx)


text:
# TODO: somehow merge the common codes below

"""
Approximate the summation below using single sample Russian roulette sampling.

S = \\sum_{i=1}^{\\infty} X_i
"""
roll(X::Function, p::DiscreteUniDist, n_mc::Int=1) = roll(X, p, rand(p, n_mc))
function roll(X::Function, p::DiscreteUniDist, τs::AbstractVector{Int})
    n_mc = length(τs)
    τs = filter(t -> t > 0, τs)
    if length(τs) == 0
        return 0
    end
    τ_min, τ_max = extrema(τs)
    ws = 1 ./ ccdf.(Ref(p), 0:τ_max-1)
    if τ_max > τ_min
        rws = reverse(cumsum(reverse(counts(τs)))) / n_mc
        ws[τ_min:τ_max] .*= rws
    end
    return sum(X.(1:τ_max) .* ws)
end

"""
Approximate the summation below using single sample Russian roulette sampling.

S = \\sum_{i=1}^{\\infty} m_i T_i
"""
roll_expectation(T::Function, m::DiscreteUniDist, p::DiscreteUniDist, n_mc::Int=1) = roll_expectation(T, m, p, rand(p, n_mc))
function roll_expectation(T::Function, m::DiscreteUniDist, p::DiscreteUniDist, τs::AbstractVector{Int})
    n_mc = length(τs)
    τs = filter(t -> t > 0, τs)
    if length(τs) == 0
        return 0
    end
    τ_min, τ_max = extrema(τs)
    ws = pdf.(Ref(m), 1:τ_max) ./ ccdf.(Ref(p), 0:τ_max-1)
    if τ_max > τ_min
        rws = reverse(cumsum(reverse(counts(τs)))) / n_mc
        ws[τ_min:τ_max] .*= rws
    end
    return sum(T.(1:τ_max) .* ws)
end

roll_expectation(T::Function, p::LogitNPD, n_mc::Int=1) = roll_expectation(T, p, rand(p, n_mc))
function roll_expectation(T::Function, p::LogitNPD, τs::AbstractVector{Int})
    n_mc = length(τs)
    τ_min, τ_max = extrema(τs)
    ws = 1 .- getρ.(Ref(p), 1:τ_max)
    if τ_max > τ_min
        rws = reverse(cumsum(reverse(counts(τs)))) / n_mc
        ws[τ_min:τ_max] .*= rws
    end
    return sum(T.(1:τ_max) .* ws)
end


text:
exp_name = "AL-kappa"

data_dirs = ["Ionosphere","PageBlocks"]

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)
oracle_param = Dict{Symbol, Any}( :type => PoolOracle, :param => Dict{Symbol, Any}())

quality_metrics = Dict(:kappa => cohens_kappa, :f1 => f1_score)
selection_metric = :kappa

n_iter = 50
n_c_steps = 20

learning_methods = vcat(
    # Kernel Active Learning with custom selection strategy
    [ALKernel(nn_stability = nns,
             stability_window = 10,
             stability_threshold = st,
             C_steps = n_c_steps,
             max_iteration = n_iter,
             initial_n = 4,
             max_gamma = 100.0,
             adjusted = true) for nns in 1:1:30 for st in [0.0]])

seeds = [i * 42 for i in 1:10]

max_stars_repo_name: holtri/lama

max_issues_repo_name: holtri/lama

max_forks_repo_name: holtri/lama