text:
module Keyboard

export readsingleline, readkey, queryrepeat, queryyesno, helpkeyboard,
        readradiomenu, readmultimenu, setmargin

import REPL
using Printf
using InteractiveUtils

# cf https://stackoverflow.com/questions/60954235/how-can-i-test-whether-stdin-has-input-available-in-julia

"""
Global var setting the left offset in the REPL for printing the list of keys/descriptors.

To modify, call setmargin().
"""
MARGIN = 1

ALLSPECIALKEYS_LUT = Dict{UInt64, Symbol}(
    0x1B5B5B41 => :FN1,
    0x1B5B5B42 => :FN2,
    0x1B5B5B43 => :FN3,
    0x1B5B5B44 => :FN4,
    0x1B5B5B45 => :FN5,
    0x1B5B31377E => :FN6,
    0x1B5B31387E => :FN7,
    0x1B5B31397E => :FN8,
    0x1B5B32307E => :FN9,
    0x1B5B32317E => :FN10,
    0x1B5B32337E => :FN11, # = SHIFT_FN1
    0x1B5B32347E => :FN12, # = SHIFT_FN2
    # 0x1B5B32337E => :SHIFT_FN1,
    # 0x1B5B32347E => :SHIFT_FN2,
    # 0x1B5B32357E => :SHIFT_FN3,
    # 0x1B5B32367E => :SHIFT_FN4,
    # 0x1B5B32387E => :SHIFT_FN5,
    # 0x1B5B32397E => :SHIFT_FN6,
    # 0x1B5B33317E => :SHIFT_FN7,
    # 0x1B5B33327E => :SHIFT_FN8,
    # 0x1B5B33337E => :SHIFT_FN9,
    # 0x1B5B33347E => :SHIFT_FN10,
    # 0x1B5B323324 => :SHIFT_FN11,
    # 0x1B5B323424 => :SHIFT_FN12,
    0x1B5B31315E => :CTRL_FN1,
    0x1B5B31325E => :CTRL_FN2,
    0x1B5B31335E => :CTRL_FN3,
    0x1B5B31345E => :CTRL_FN4,
    0x1B5B31355E => :CTRL_FN5,
    0x1B5B31375E => :CTRL_FN6,
    0x1B5B31385E => :CTRL_FN7,
    0x1B5B31395E => :CTRL_FN8,
    0x1B5B32305E => :CTRL_FN9,
    0x1B5B32315E => :CTRL_FN10,
    0x1B5B32335E => :CTRL_FN11,
    0x1B5B32345E => :CTRL_FN12,
    0x1B1B5B5B41 => :ALT_FN1,
    0x1B1B5B5B42 => :ALT_FN2,
    0x1B1B5B5B43 => :ALT_FN3,
    0x1B1B5B5B44 => :ALT_FN4,
    0x1B1B5B5B45 => :ALT_FN5,
    0x1B1B5B31377E => :ALT_FN6,
    0x1B1B5B31387E => :ALT_FN7,
    0x1B1B5B31397E => :ALT_FN8,
    0x1B1B5B32307E => :ALT_FN9,
    0x1B1B5B32317E => :ALT_FN10,
    0x1B1B5B32337E => :ALT_FN11,
    0x1B1B5B32347E => :ALT_FN12,
    0x1B5B337E => :DEL,
    0x1B5B317E => :HOME,
    0x1B5B347E => :END,
    0x1B5B357E => :PAGEUP,
    0x1B5B367E => :PAGEDOWN,
    0x1B5B333B327E => :SHIFT_DEL,
    0x1B5B313B327E => :SHIFT_HOME,
    0x1B5B343B327E => :SHIFT_END,
    0x1B5B353B327E => :SHIFT_PAGEUP,
    0x1B5B363B327E => :SHIFT_PAGEDOWN,
    0x1B5B333B357E => :CTRL_DEL,
    0x1B5B313B357E => :CTRL_HOME,
    0x1B5B343B357E => :CTRL_END,
    0x1B5B353B357E => :CTRL_PAGEUP,
    0x1B5B363B357E => :CTRL_PAGEDOWN,
    0x1B1B5B337E => :ALT_DEL,
    0x1B1B5B317E => :ALT_HOME,
    0x1B1B5B347E => :ALT_END,
    0x1B1B5B357E => :ALT_PAGEUP,
    0x1B1B5B367E => :ALT_PAGEDOWN,
    0x1B5B41 => :CURSORUP,
    0x1B5B42 => :CURSORDOWN,
    0x1B5B44 => :CURSORLEFT,
    0x1B5B43 => :CURSORRIGHT,
    0x1B5B313B3541 => :CTRL_CURSORUP,
    0x1B5B313B3542 => :CTRL_CURSORDOWN,
    0x1B5B313B3544 => :CTRL_CURSORLEFT,
    0x1B5B313B3543 => :CTRL_CURSORRIGHT,
    0x1B1B5B41 => :ALT_CURSORUP,
    0x1B1B5B42 => :ALT_CURSORDOWN,
    0x1B1B5B44 => :ALT_CURSORLEFT,
    0x1B1B5B43 => :ALT_CURSORRIGHT
    )

CTRLCHARS_LUT =  Dict(zip(collect(Char(1):Char(31)), collect("ABCDEFGHIJKLMNOPQRSTUVWXYZ[\\]^_")))

"""
            _vector_sans_type_tostring(vec::AbstractVector)::String

Return a String detailing the arg 'vec''s values with nesting (if any), without the Types moiety.
"""
function _vector_sans_type_tostring(vec::AbstractVector)::String
    io = IOBuffer()

    function recurse(v::AbstractVector)
        print(io, "[")
        for (i, elt) in enumerate(v)
            i > 1 && print(io, ", ")
            if elt isa AbstractVector
                recurse(elt)
            else
                if elt isa Symbol
                    elt = String(":" * String(elt))
                elseif elt isa Char
                    elt = "'" * elt * "'"
                elseif elt isa String
                    elt = "\"" * elt * "\""
                end
                print(io, elt)
            end
        end
        print(io, "]")
    end
    recurse(vec)

    return String(take!(io))
end

"""
            _splitspecialkeycombo(s::Symbol) ::
                    NamedTuple{(:metakey, :specialkey), Tuple{Union{Symbol, Nothing}, Symbol}}

Split the supplied Symbol, arg 's' according to the underscore seperator '_'.
The left-hand moiety will be a metakey name; the right a specialkey name. If no seperator is present, assume 's'
is a specialkey.

# Examples of arg 's':

 :CTRL_HOME, :HOME

 Returns (metakey=:CTRL, specialkey=:HOME) and (metakey=nothing, specialkey=:HOME)
"""
function _splitspecialkeycombo(s::Symbol) ::
        NamedTuple{(:metakey, :specialkey), Tuple{Union{Symbol, Nothing}, Symbol}}
    ans1 = nothing
    ans2 = nothing
    a = split(String(s), "_")
    if length(a) == 1
        ans2 = Symbol(a[1])
    else
        ans1 = Symbol(a[1])
        ans2 = Symbol(a[2])
    end

    return (metakey = ans1, specialkey = ans2)
end

"""
            _getallkeynames()::
                    Tuple{NamedTuple{(:ctrl, :shift, :alt), Tuple{Symbol, Symbol, Symbol}}, NTuple{21,Symbol}}
# Return
All extracted meta- and special-keys as symbols.

(The meta-keys are in the NamedTuple, eg alt=:ALT, and the special-keys are in the NTuple, eg :ESC).
"""
function _getallkeynames() ::
        Tuple{NamedTuple{(:ctrl, :shift, :alt), Tuple{Symbol, Symbol, Symbol}}, NTuple{21,Symbol}}
    metas = Set{Symbol}()
    specials = Set{Symbol}()
    for v in values(ALLSPECIALKEYS_LUT)
        a = _splitspecialkeycombo(v)
        push!(specials, a.specialkey)
        if a.metakey !== nothing
            push!(metas, a.metakey)
        end
    end
    if length(metas) != 3
        error("Metakey-name error in ALLSPECIALKEYS_TABLE")
    end
    if length(specials) != 21
        error("Special-name error in ALLSPECIALKEYS_TABLE")
    end
    metakeys = collect(metas)
    sort!(metakeys)
    ans1 = (ctrl = metakeys[2], shift = metakeys[3], alt = metakeys[1])
    specialsnames = collect(specials)
    sort!(specialsnames)

    return ans1, tuple(specialsnames ...)
end

"""
            _isprintable(c::Char) :: Bool = c > ' ' && isprint(c)
Here, printable means that the arg 'c' > space and isprint(c) evaluates to true.
"""
_isprintable(c::Char) :: Bool = c > ' ' && isprint(c)
"""
            _isprintable(s::AbstractString) :: Bool
Test arg 's' consists of _isprintable() characters (space is also regarded as non-printable).
"""
function _isprintable(s::AbstractString) :: Bool
    isempty(s) && return false
    for c in s
        !_isprintable(c) && return false
    end
    return true;
end

METAKEYNAMES_TABLE, SPECIALKEYNAMES_TABLE = _getallkeynames()

CTRLSPECIALKEYNAMES_LUT = Dict{Char, Symbol}(
    '[' => :ESC,
    'H' => :BACK,
    'I' => :TAB,
    'M' => :ENTER,
)

"""
            readkey(; prompt::String = "Press any key..")

If arg 'prompt' != "" then the function evaluates println(prompt) before waiting for a key-press.

        Returns a NamedTuple{(:metakey, :keychar, :specialkey),
            Tuple{Union{Symbol, Nothing}, Union{Char, Nothing}, Union{Symbol, Nothing}}}
"""
function readkey(; prompt::String="Press any key..") ::
    NamedTuple{(:metakey, :keychar, :specialkey),
        Tuple{Union{Symbol, Nothing}, Union{Char, Nothing}, Union{Symbol, Nothing}}}
    term = REPL.Terminals.TTYTerminal("xterm",stdin,stdout,stderr)
    REPL.Terminals.raw!(term,true)
    Base.start_reading(stdin)

    prompt != "" && println(prompt)
    read(stdin, bytesavailable(stdin)) # clear stdin

    while true
        keychar::Union{Char, Nothing} = nothing
        metakey::Union{Symbol, Nothing} = nothing
        specialkey::Union{Symbol, Nothing} = nothing

        sleep(0.01)
        ba = bytesavailable(stdin)
        ba == 0 && continue

        data = read(stdin, ba)
        if ba == 1
            keychar = data[1]
            if iscntrl(keychar)
                keychar = CTRLCHARS_LUT[keychar]
                metakey = METAKEYNAMES_TABLE.ctrl
                if haskey(CTRLSPECIALKEYNAMES_LUT, keychar)
                    specialkey = CTRLSPECIALKEYNAMES_LUT[keychar]
                end
            end
        elseif ba == 2
            keychar = data[2]
            metakey = METAKEYNAMES_TABLE.alt
        else
            n::UInt64 = data[1]
            for i in 2:ba
                n <<= 8
                n += data[i]
            end
            if haskey(ALLSPECIALKEYS_LUT, n)
                combo = ALLSPECIALKEYS_LUT[n]
                metakey, specialkey = _splitspecialkeycombo(combo)
                keychar = nothing
            else
                continue # key not in LUTs
            end
        end
        metakey === keychar === specialkey === nothing && continue
        ans = (metakey = metakey, keychar = keychar, specialkey = specialkey)

        return ans
    end
end

"""
Compose and return a user-friendly string describing special key choices to be printed after the list of choices upon return.
# Examples
        "[ESC=cancel, Enter=y, ^C=quit]" :- escape = quitCtrlC = true; defaultKey isa Char.
        "[ESC=cancel, ^C=quit, TAB=all, BACK=none, ENTER=commit]" :- escape = quitCtrlC = ismultisel = true.
"""
function _getescenter(escape::Bool, defaultKey::Union{Char, Nothing}, quitCtrlC::Bool, ismultisel::Bool)::String
    s = ""
    function escenterquit()
        escape && (s = "[ESC=cancel")
        if defaultKey === nothing && !quitCtrlC
            if isempty(s)
                return s
            else
                return s * ']'
            end
        end
        s = isempty(s) ? "[" : s
        if defaultKey isa Char
            if escape
                s *= ", Enter=$defaultKey"
            else
                s *= "Enter=$defaultKey"
            end
        end
        !quitCtrlC && return s * ']'
        if length(s) > 2
            s *= ", ^C=quit]"
        else
            s = "[^C=quit]"
        end
        return s
    end

    s = escenterquit()
    !ismultisel && return s

    s = strip(s) # just in case
    s = isempty(s) ? "[" : s[1:end-1] * ", "
    s *= "TAB=all, BACK=none, ENTER=commit]"
    return s
end

"""
            _getkeyboardchar(
                keys::String, # a string of Chars with which to validate a key-press.
                defaultKey::Union{Char, Nothing}, # if this is a Char, detect the Enter key press.
                case::Bool, # whether to apply case-sensitivity to detection in arg 'keys'.
                escape::Bool, # whether to detect Esc key-press.
                quitCtrlC::Bool, # whether to detect Ctrl+C key-press.
                ismultisel::Bool # true = detect Tab (fill), Back (unfill) and Enter (commit) key-presses.
                )::Tuple{Union{Char, Symbol}, Int64}

Wait for a key-press. Validate using the args supplied.
# Return a Tuple{Union{Char, Symbol}, Int64} upon valid key-press:
        1) if key-press = Ctrl+C, return (Char(3), index = -1).
        2) if key-press is in arg 'keys' then return (keychar, index_into_keys).
        3) if key-press is a special key then return (specialkey, index = 0).
"""
function _getkeyboardchar(
    keys::String, # a list of Chars with which to validate a key-press
    defaultKey::Union{Char, Nothing}, # if this is a Char, detect the Enter key press
    case::Bool, # whether to apply case-sensitivity to detection in arg 'keys'
    escape::Bool, # whether to detect Esc key-press
    quitCtrlC::Bool, # whether to detect Ctrl+C key-press
    ismultisel::Bool # true = detect Tab (fill), Back (unfill) and Enter (commit) key-presses
    ) :: Tuple{Union{Char, Symbol}, Int64}

    ix = nothing
    u = nothing
    while true
        res = readkey(prompt="")
        quitCtrlC && res.metakey === METAKEYNAMES_TABLE.ctrl && res.keychar == 'C' &&
                (u=Char(3); ix=-1; break)
        escape && res.specialkey === :ESC &&
                (u=:ESC; ix=0; break)
        defaultKey !== nothing && res.specialkey === :ENTER &&
                (u=defaultKey; ix=findfirst(u, keys); break)
        sk = res.specialkey
        ismultisel && sk in (:TAB, :ENTER, :BACK) &&
                (u = sk; ix=0; break)
        (res.keychar === nothing || res.metakey !== nothing) && continue
        u = res.keychar
        !case && (u = lowercase(u))
        ix = findfirst(u, keys)
        ix !== nothing && break
    end
    return u, ix
end

"""
        readsingleline(
           prompt::String,
           keys::String, # a list of keyboard characters one of which the user may select (not case-sensitive).
           ;
           escape::Bool = true,
           defaultKey::Union{Char, Nothing} = nothing, # or a keyboard character present in the arg 'keys' (not case-sensitive).
           quitCtrlC::Bool = false
           )::Union{Char, Nothing}

Prints the arg 'prompt' followed by waiting for a valid key-press.

It returns EITHER a user-selected char (always in lower case) present in the arg 'keys' OR nothing (the user pressed ESC).

If the arg 'quitCtrlC' = true THEN returns Char(3) if Ctrl+C is pressed.

NB. This function is not case-sensitive.
"""
function readsingleline(
    prompt::String,
    keys::String,
    ;
    escape::Bool = true,
    defaultKey::Union{Char, Nothing} = nothing,
    quitCtrlC::Bool = false
    )::Union{Char, Nothing}

    keychars = Vector{Char}("(")
    prompt = strip(prompt)
    isempty(prompt) && error("The arg 'prompt' has no text.")
    keys = lowercase(join(split(keys)))
    isempty(keys) && error("The arg 'keys' is empty.")
    !_isprintable(keys) && error("The arg 'keys' (\"$keys\") must consist of printable characters.")
    if defaultKey isa Char
        defaultKey = lowercase(defaultKey)
        !(defaultKey in keys) && error("The arg 'defaultKey'('$defaultKey') is not in the arg 'keys'(\"$keys\").")
    end
    for chr in keys
        push!(keychars, chr, '|')
    end
    pop!(keychars); push!(keychars, ')')
    prompt *= join(keychars) * _getescenter(escape, defaultKey, quitCtrlC, false)

    println(prompt)

    chr, _ = _getkeyboardchar(keys, defaultKey, false, escape, quitCtrlC, false)
    chr === :ESC && return nothing
    return chr
end

"""
        _readmultiline(
            prompt::String,  # a string to print first of all, describing the menu.
            data,  # a Vector of Vectors with which to extract keys (optional), descriptors and functions (optional).
            ;
            escape::Bool=true,
            case::Bool=false, # true=case-sensitive for the extracted keys.
            ismultisel::Bool, # true=multiple keys selection; else radio type selection
            quitCtrlC::Bool=false
            )::Union{T, Union{Tuple{Vararg{T}}, Tuple{}}, Nothing, Char}
                    where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}

# Argument 'data'

            Examples:

            1) The descriptors are mandatory; keys will be auto-generated:
                [
                ["aaa"],
                ...,
                ["eee"],
                ]
            2) Keys + descriptors:
                [
                ['z',"aaa"],
                ...,
                ['u',"eee"],
                ]
            3) Keys + descriptors + functions (zero-arg):
                [
                ['z',"aaa", fun_1],
                ...,
                ['u',"eee", fun_5],
                ]
            4) Descriptors + functions (zero-arg):
                [
                ["aaa", fun_1],
                ...,
                ["eee", fun_5],
                ]
            Note that
                (1) dummy()=nothing  may be used as a dummy function.
                (2) If functions are specified in a list, all items must contain one.

# 4 types of return:

        1) a Char = Char(3), Ctrl+c was pressed.
        2) a Nothing - ESC key was pressed.
        3) a NamedTuple{(:keychar, :index), Tuple{Char, Int64}} - a key character was selected;
            :index is the position in the list of keys (ismultisel=false).
        4) a Union{Tuple{Vararg{T}}, Tuple{}} where T <:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}
            This is either a tuple of T's above, or an empty list; (ismultisel=true).
"""
function _readmultiline(
    prompt::String,
    data,
    ;
    escape::Bool=true,
    case::Bool=false,
    ismultisel::Bool,
    quitCtrlC::Bool=false
    )::Union{T, Union{Tuple{Vararg{T}}, Tuple{}}, Nothing, Char} where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}

    datauniform = nothing
    try; datauniform = Vector{Vector{Any}}(deepcopy(data))
    catch; error("The passed 'data' must be of the form [[...],[...], etc]")
    end

    sel = ("$(' '^MARGIN){", "} ") # prefix the keychar, and postfix for a user-selected key
    unsel = ("$(' '^MARGIN) ", "  ") # ditto for a non-user-selected key

    # ASSIGN vars typeel_1 AND typeel_3 (CORRESPONDING TO THE TYPES OF THE 3 COLUMNS IN datauniform)
    # AND ASSIGN nothing TO UNASIGNED COLUMNS IN datauniform:
    typeel_1 = Nothing; typeel_3 = Nothing
    firstitem = data[1]
    lenfirstitem = length(firstitem)

    if lenfirstitem == 0
        error("The first item is empty in the passed arg 'data'")
    elseif lenfirstitem == 1
        push!.(datauniform, nothing)
        pushfirst!.(datauniform, nothing)
    elseif lenfirstitem == 2
        fun = pushfirst!
        if firstitem[2] isa String
            fun = push!
            typeel_1 = Char
        else
            typeel_3 = Function
        end
        fun.(datauniform, nothing)
    elseif lenfirstitem == 3
        typeel_1 = Char
        typeel_3 = Function
    else
        error("The number of elements can not exceed 3 in the item: $(_vector_sans_type_tostring(firstitem))")
    end

    typeel_1 === Nothing && (case = false) # auto-generated keys are all lower case anyway

    # EXTRACT ALL ELEMENTS FROM datauniform INTO 3 Vectors:
    keychars = Vector{typeel_1}()
    descriptors = Vector{String}()
    functions = Vector{typeel_3}()
    for i in 1:length(datauniform)
        item = datauniform[i]
        dataitem = data[i]
        length(dataitem) != lenfirstitem &&
            error("There must be $lenfirstitem element(s) in $(_vector_sans_type_tostring(dataitem))")
        err = nothing
        try
            err = "The first element must be a Char in "
            ele = item[1] isa typeel_1 ? item[1] : "dummy" # int and floats are converted to Char in push!
            ele isa Char && !_isprintable(ele) &&
                (err = "The key character must be printable in "; error())
            push!(keychars, ele)
            err = "Either a String is missing or it is in the wrong position in "
            push!(descriptors, item[2])
            err = "The last element must be a Function in "
            push!(functions, item[3])
        catch
            error(err * _vector_sans_type_tostring(dataitem))
        end
    end

    # CHECK FOR DUPLICATES IN keychars AND CONVERT IT TO A CASE-ADJUSTED String keystring:
    keystring = nothing
    if typeel_1 === Char
        keystring = join(keychars)
        !case && (keystring = lowercase(keystring))
        length(Set(collect(keystring))) < length(keystring) &&
            error("There are duplicates in the extracted keys: \"$(join(keychars))\"")
    else
        keystring = join('1':'9') * join('a':'z')
        length(data) > length(keystring) &&
            error("Can't auto-generate all of the key characters - the arg 'data' is too big.")
        keystring = keystring[1:length(data)]
    end

    # OK, NO DETECTED ERRORS IN PASSED 'data'; PRINT THE LIST FOLLOWED BY THE PROMPT:
    println(prompt)
    for i in 1:length(data)
        println("$(unsel[1])$(keystring[i])$(unsel[2])$(descriptors[i])")
    end
    println(_getescenter(escape, nothing, quitCtrlC, ismultisel))

    # DETECT KEY-PRESSES AND ASSIMILATE THEM:
    lenkeys = length(keystring)
    selkeys = Vector{Bool}(undef, lenkeys)
    fill!(selkeys, false)
    while true
        cors, index = _getkeyboardchar(keystring, nothing, case, escape, quitCtrlC, ismultisel)
        cors === Char(3) && return cors
        cors === :ESC && return nothing
        if !ismultisel
            typeel_3 === Function && cors isa Char && cors !== Char(3) && functions[index]() # eval the function
            return (keychar=cors, index=index)
        end
        if cors === :ENTER
            v = Vector{NamedTuple{(:keychar, :index), Tuple{Char, Int64}}}()
            for i in 1:lenkeys
                selkeys[i] && push!(v, (keychar=keystring[i], index=i))
                typeel_3 === Function && selkeys[i] && functions[i]() # eval a function
            end
            return Tuple(v)
        elseif cors === :TAB
            fill!(selkeys, true)
        elseif cors === :BACK
            fill!(selkeys, false)
        else
            selkeys[index] = !selkeys[index] # keychar pressed - toggle selection
        end
        print("\x1b[?25l\x1b[$(lenkeys + 1)A") # move cursor up n lines and hide it
        for i in 1:lenkeys
            p = selkeys[i] ? sel : unsel
            println("$(p[1])$(keystring[i])$(p[2])")
        end
        println("\x1b[?25h") # extra 1 for the keys info line at the bottom; show cursor
    end
end

"""
            readradiomenu(
                    prompt::String, # a string to print first of all, describing the menu.
                    data::Vector{String}, # a single Vector of Strings
                    ;
                    escape::Bool=true,
                    case::Bool=false, # true=case-sensitive for the extracted keys.
                    quitCtrlC::Bool=false
                    )::Union{T, Nothing, Char} where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}

This overloaded function differs in that the arg 'data' must be a Vector of Strings.
"""
function readradiomenu(
        prompt::String, # a string to print first of all, describing the menu.
        data::Vector{String}, # a single array of Strings
        ;
        escape::Bool=true,
        case::Bool=false, # true=case-sensitive for the extracted keys.
        quitCtrlC::Bool=false
        )::Union{T, Nothing, Char} where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}

    d = [[x] for x in data]
    readradiomenu(
        prompt,
        d,
        ;
        escape = escape,
        case = case,
        quitCtrlC = quitCtrlC
        )
end

"""
            readradiomenu(
                prompt::String, # a string to print first of all, describing the menu.
                data, # a Vector of Vectors with which to extract keys (optional), descriptors and functions (optional).
                ;
                escape::Bool=true,
                case::Bool=false, # true=case-sensitive for the extracted keys.
                quitCtrlC::Bool=false
                )::Union{T, Nothing, Char} where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}

Prints a prompt followed by a formatted list of keys and descriptors; it then waits for a valid single keypress,
upon which the function will return with the results. If function names were supplied, the corresponding zero-arg
function will be evaluated.

Evaluate setmargin(..) to set the left offset into the list printout (global).

# Argument 'data'

        Examples:

        1a) The descriptors are mandatory; keys will be auto-generated:
            ["aaa", ..., "eee"]
        1b) The descriptors are mandatory; keys will be auto-generated:
            [
            ["aaa"],
            ...,
            ["eee"],
            ]
        2) Keys + descriptors:
            [
            ['z',"aaa"],
            ...,
            ['u',"eee"],
            ]
        3) Keys + descriptors + functions (zero-arg):
            [
            ['z',"aaa", fun_1],
            ...,
            ['u',"eee", fun_5],
            ]
        4) Descriptors + functions (zero-arg); keys will be auto-generated:
            [
            ["aaa", fun_1],
            ...,
            ["eee", fun_5],
            ]
        Note that
            (1) dummy()=nothing  may be used as a dummy function.
            (2) If functions are specified in a list, all items must contain one.

# 3 types of return:

    1) a Char = Char(3) - Ctrl+C was pressed.
    2) a Nothing - the ESC key was pressed.
    3) a NamedTuple{(:keychar, :index), Tuple{Char, Int64}} - a key character was selected;
                                :index is the position in the list of keys/descriptors.
"""
readradiomenu(
    prompt::String,
    data,
    ;
    escape::Bool=true,
    case::Bool=false,
    quitCtrlC::Bool=false
    )::Union{T, Nothing, Char} where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}} =
        _readmultiline(prompt, data, escape=escape, case=case, ismultisel=false, quitCtrlC=quitCtrlC)

"""
            readmultimenu(
                prompt::String, # a string to print first of all, describing the menu.
                data::Vector{String}, # a single Vector of Strings
                ;
                escape::Bool=true,
                case::Bool=false, # true=case-sensitive for the extracted keys.
                quitCtrlC::Bool=false
                )::Union{Union{Tuple{Vararg{T}}, Tuple{}}, Nothing, Char}
                        where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}

This overloaded function differs in that the arg 'data' must be a Vector of Strings.
"""
function readmultimenu(
    prompt::String, # a string to print first of all, describing the menu.
    data::Vector{String}, # a single array of Strings
    ;
    escape::Bool=true,
    case::Bool=false, # true=case-sensitive for the extracted keys.
    quitCtrlC::Bool=false
    )::Union{Union{Tuple{Vararg{T}}, Tuple{}}, Nothing, Char} where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}

    d = [[x] for x in data]
    readmultimenu(
        prompt,
        d,
        ;
        escape = escape,
        case = case,
        quitCtrlC = quitCtrlC
        )
end

"""
            readmultimenu(
                prompt::String,  # a string to print first of all, describing the menu.
                data,  # a Vector of Vectors with which to extract keys (optional), descriptors and functions (optional).
                ;
                escape::Bool=true,
                case::Bool=false, # true=case-sensitive for the extracted keys.
                quitCtrlC::Bool=false
                )::Union{Union{Tuple{Vararg{T}}, Tuple{}}, Nothing, Char}
                    where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}

Prints a prompt followed by a formatted list of keys and descriptors;
it then waits for a succession of valid single keypresses, upon which each selected key from the list
will be 'un|checked' (toggled). If function names were supplied, the corresponding zero-arg functions will
be evaluated sequentially. Pressing the Enter key will commit the changes.

Pressing TAB will check all of the keys; whilst pressing BACK will uncheck them.

Evaluate setmargin(..) to set the left offset into the list's printout (global).

# Argument 'data'

            Examples:

            1a) The descriptors are mandatory; keys will be auto-generated:
                [ "aaa", ..., "eee"]
            1b) The descriptors are mandatory; keys will be auto-generated:
                [
                ["aaa"],
                ...,
                ["eee"],
                ]
            2) Keys + descriptors:
                [
                ['z',"aaa"],
                ...,
                ['u',"eee"],
                ]
            3) Keys + descriptors + functions (zero-arg):
                [
                ['z',"aaa", fun_1],
                ...,
                ['u',"eee", fun_5],
                ]
            4) Descriptors + functions (zero-arg); keys will be auto-generated:
                [
                ["aaa", fun_1],
                ...,
                ["eee", fun_5],
                ]
            Note that
                (1) dummy()=nothing  may be used as a dummy function.
                (2) If functions are specified in a list, all items must contain one.

# 3 types of return:

        1) a Char = Char(3), Ctrl+c was pressed.
        2) a Nothing - ESC key was pressed.
        3) a Union{Tuple{Vararg{T}}, Tuple{}} where T <:NamedTuple{(:keychar, :index), Tuple{Char, Int64}}
            This is either a tuple of T's above, or an empty tuple.
"""
readmultimenu(
    prompt::String,
    data,
    ;
    escape::Bool=true,
    case::Bool=false,
    quitCtrlC::Bool=false
    )::Union{Union{Tuple{Vararg{T}}, Tuple{}}, Nothing, Char} where T>:NamedTuple{(:keychar, :index), Tuple{Char, Int64}} =
        _readmultiline(prompt, data, escape=escape, case=case, ismultisel=true, quitCtrlC=quitCtrlC)

"""
                queryrepeat(prompt::String="Repeat?") :: Bool

Print the prompt and wait for a valid key press; upon which the function will return a Bool.

IF return = true, the 'y' key or Enter key was pressed ELSE the 'n' key or ESC key was pressed.
"""
function queryrepeat(prompt::String="Repeat?") :: Bool
    c = readsingleline(prompt, "YN", escape=true, defaultKey='Y')
    (c === nothing || c === 'n') && return false
    return true
end

"""
        queryyesno(prompt::String)::Bool

Print the prompt and wait for a valid key press; upon which the function will return a Bool.

IF return = true, the 'y' key or Enter key was pressed ELSE the 'n' key or ESC key was pressed.
"""
queryyesno(prompt::String)::Bool = queryrepeat(prompt)

"""
            helpkeyboard()

Print all of the special key names extracted (including the metakeys).
"""
function helpkeyboard()
    f(x) = ":$(String(x))"
    L = f.(SPECIALKEYNAMES_TABLE) # convert tuple of symbols to tuple of strings
    L2 = [s for s in L]
    sort!(L2)
    functionkeys = [s for s in L2 if lowercase(s[2]) == 'f']
    specialkeys = [s for s in L2 if lowercase(s[2]) != 'f']
    metakeys = f.(Tuple(METAKEYNAMES_TABLE))
    ctrlspecialkeys = f.(values(CTRLSPECIALKEYNAMES_LUT))
    step = 6
    function printit(title, list)
        println(title)
        for i in 0:step:length(list) - 1
            for j in 1:step
                ix = i + j
                ix > length(list) && break
                print("\t$(list[ix])")
            end
            println()
        end
    end
    printit("\nFunction-key symbols:", functionkeys)
    printit("\nSpecial-key symbols:", specialkeys)
    printit("\nMeta-key symbols:", metakeys)
    printit("\nControl-special-key symbols:", ctrlspecialkeys)
end

"""
            setmargin(margin::Real)

Sets the private global MARGIN which is used to set the left offset of the printed list of keys/descriptors.
"""
function setmargin(margin::Real)
    x::Int64 = floor(Int64, margin)
    global MARGIN = min(max(x, 0), 30)
end

end # Keyboard


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule timg_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("timg")
JLLWrappers.@generate_main_file("timg", UUID("6c2e2082-7334-5db5-9567-78a40bc96671"))
end  # module timg_jll

ext: jl


text:
using QF939
using Test

@testset "Distance function" begin

	x = [0, 0];
	y = [12, 12];
	side = 10;

	@test isapprox(QF939.rpbc(x, y, side), sqrt(8));

	y = [17, 17];

	@test isapprox(QF939.rpbc(x, y, side), sqrt(18));
end

@testset "Potential Energy" begin

	data = Data();

	@test isapprox(QF939.upair(π, data), -0.0003250451775148634);

	p = initial_point(data);

	@test isapprox(utotal(p, data), -7472.164849146468);
end

@testset "Potential Energy and Forces" begin

	data = Data();
	p = initial_point(data);

	u, f = uftotal(p, data);

	@test isapprox(u, -7472.164849146468);
	@test isapprox(f[1], [2.9749744219446326, 1.730790463586403]);
	@test isapprox(f[data.N], [-3.1208157043186247, -2.0513425336203057]);
end


text:
module PMachine

export Machine, exe, load

# N steht fUr "numerischer" Typ,
# also integer, real oder Adresse, T fur beliebigen einfachen Typ, also fur numerische,
# logische, Zeichen- und Aufzahlungstypen und fur Adressen, i fiir integer, r
# fUr real, b fUr boolean, a fur Adresse

# R = Float64
# B = Bool
# A = Int64
# I = Int64
#
# N = Union{ I, R, A }
# T = Union{ N, B }
#
# C = Union{ Symbol, T }
#
# mutable struct Machine
#     SP      :: A
#     store   :: Vector{ T }
#     PC      :: A
#     code    :: Vector
# end
# Machine( code::Vector = [] ) = Machine( 0, zeros( 100 ), 0, code )
#
# function reset( m::Machine )
#     m.code = []
#     m.PC = 0
# end
#
# function load( m::Machine, c::Vector )
#     reset( m )
#     m.code = c;
# end
#
# function next_code( m::Machine )
#     m.PC += 1
#     m.code[ m.PC ]
# end
#
# function exe_type_symbol_num( m::Machine )
#     t = next_code( m )
#     if typeof( t ) != Symbol
#         throw( "No type symbol after :add, $(t)" )
#     end
#     if t == :i
#         ty = I
#     elseif t == :r
#         ty = R
#     elseif t == :a
#         ty = A
#     else
#         throw( "Unknown numeric type symbol, $( t )" )
#     end
#     ty
# end
#
# function exe_add( m::Machine )
#     ty =  exe_type_symbol_num( m )
#     left  = m.store[ m.SP - 1 ]
#     right = m.store[ m.SP ]
#     if typeof( left ) != ty
#         throw( "Wrong type on SP - 1  for :add, $( left )")
#     end
#     if typeof( right ) != ty
#         throw( "Wrong type on SP  for :add, $( right )")
#     end
#     println( "ADD")
#     m.store[ m.SP - 1 ] = left + right
#     m.SP -= 1
# end
#
# function exe_type_symbol( m::Machine )
#     t = next_code( m )
#     if typeof( t ) != Symbol
#         throw( "No type symbol after :add, $(t)" )
#     end
#     if t == :i
#         ty = I
#     elseif t == :r
#         ty = R
#     elseif t == :a
#         ty = A
#     elseif t == :b
#         ty = B
#     else
#         throw( "Unknown type symbol, $( t )" )
#     end
#     ty
# end
#
#
# function exe_ldc( m::Machine )
#     ty =  exe_type_symbol( m )
#     imm = next_code( m )
#     if typeof( imm ) != ty
#         throw( "Wrong type for :ldc, $( imm )")
#     end
#     m.SP += 1
#     m.store[ m.SP ] = imm
# end
#
#
#
# function exe( m::Machine )
#     run = true
#     while run
#         s = next_code( m )
#         println( s )
#         if s == :nop
#             println( "NOP!" )
#         elseif s == :add
#             exe_add( m )
#         elseif s == :ldc
#             exe_ldc( m )
#         elseif s == :halt
#             println( "SP: $( m.SP ), PC: $( m.PC )" )
#             println( "Store: " )
#             show( m.store )
#             run = false
#         end
#         # println( typeof( s ) )
#     end
# end
#
# function exe( m::Machine, c::Vector )
#     load( m, c )
#     exe( m )
# end

export Register, get, set, reset
import Base.get

mutable struct Register{ Val <: Integer } <: Integer
    value  :: Val
end
Register{T}()           where T = Register{T}( 0 )
Register( val :: T )    where T = Register{T}( val )

get( r :: Register{T} )         where T = r.value
set( r :: Register{T}, v :: T ) where T = r.value = v

reset( r :: Register{T} ) where T = set( r, T( 0 ) )

export Memory, getAddr, setAddr, read, write, reset, size

struct Memory{ Val, Add <: Integer }
    values  :: Vector{ Val }
    addr    :: Register{ Add }
end
Memory{ V, A }( size :: Integer ) where A where V = Memory{ V, A }( zeros( V, size ), Register{ A }() )
# Memory{ V, A }()                  where A where V = Memory{ V, A }( zeros( V, size ), Register{ A }() )

getAddr( m :: Memory{ V, A } )          where A where V = get( m.addr )
setAddr( m :: Memory{ V, A }, a :: A )  where A where V = set( m.addr, a )

read(  m :: Memory{ V, A } )                where A where V = m.values[ get( m.addr ) + 1 ]
write( m :: Memory{ V, A }, value :: V )    where A where V = m.values[ get( m.addr ) + 1 ] = value

size( m :: Memory{ V, A } ) where A where V = length( m.values )

function reset(  m :: Memory{ V, A } ) where A where V
    reset( m.addr )
    fill!( m.values, V( 0 ) )
    m
end

export Machine, codesize, storesize

struct Machine{ Word <: Integer }
    PC      :: Register{ Word }
    SP      :: Register{ Word }

    code    :: Memory{ Word, Word }
    store   :: Memory{ Word, Word }
end

function Machine{ W }( codesize, storesize :: Integer ) where W
    Machine{ W }(
        Register{ W }(),
        Register{ W }(),
        Memory{ W, W }( codesize ),
        Memory{ W, W }( storesize )
    )
end
Machine( codesize, storesize :: Integer ) = Machine{ Int }( codesize, storesize :: Integer )
Machine( size ) = Machine( size, size )
Machine() = Machine( 100 )

# function reset( m :: Machine{ W } ) where W
#     reset( m.PC )
#     reset( m.SP )
#     reset( m.code )
#     reset( m.store )
# end

codesize(  m :: Machine{ W } ) where W = size( m.code )
storesize( m :: Machine{ W } ) where W = size( m.store )

function load( m::Machine{ W }, c::Vector{ W }, start::W = 0 ) where W
    reset( m )
    m.code.values[ 1:end ] = c[ 1:end ];
end

# function next_code( m::Machine )
#     m.PC += 1
#     m.code[ m.PC ]
# end

end

getAddr( m :: Memory{ V, A } )          where A where V = get( m.addr )
setAddr( m :: Memory{ V, A }, a :: A )  where A where V = set( m.addr, a )

read(  m :: Memory{ V, A } )                where A where V = m.values[ getAddr( m ) + 1 ]
write( m :: Memory{ V, A }, value :: V )    where A where V = m.values[ getAddr( m ) + 1 ] = value

end


text:
function _treeplot(tree::Tree)
  nodequeue = findroots(tree)
  treesize = leafcount.(Ref(tree), nodequeue)
  distances = distance.(Ref(tree), nodequeue)
  countpct = treesize / sum(treesize)
  left_pos = cumsum(countpct) .- (0.5 * countpct)
  queueposition = 1
  while(queueposition <= length(nodequeue))
    @debug "Processing Node $(nodequeue[queueposition])"
    children = childnodes(tree, nodequeue[queueposition])
    @debug "Children Nodes = $children"
    if length(children) > 0
      append!(nodequeue, children)

      # find the proportion of all leaf nodes contained the subtrees
      subtreesize = max.(leafcount.(Ref(tree), children), 1)
      @debug "Number of leaves per child = $subtreesize"

      append!(countpct, subtreesize / sum(treesize))
      @debug "Proportion of total leaves per child = $(round.(subtreesize / sum(treesize), digits=3))"
      # calculate the positions
      # 1a - find parent node postion
      # 1b - find left-most leaf position of subtree
      subtree_left = left_pos[queueposition] - 0.5 * countpct[queueposition]
      @debug "Position of parent node = $(round(left_pos[queueposition], digits=3))"
      @debug "Leftmost descendant leaf node position = $(round(subtree_left, digits=3))"

      # 2a - find indices of children
      c_ind = (length(countpct) - length(children) + 1):length(countpct)
      @debug "Indices of children = $c_ind"
      # 2b - find position relative to left of subtree
      # 2c - use a cumulative sum based on leaf node count %
      # 2d - adjust these cumulative-sum-based positions down to their midpoints
      append!(left_pos, subtree_left .+ cumsum(countpct[c_ind]) .- 0.5 * countpct[c_ind])
      @debug "Position of children = $(round.(left_pos[c_ind], digits=3))"
    end
    queueposition += 1
    @debug "Increment queue position to $queueposition"
  end

  # distance from root + tree.height as node height/distance
  heights = distance.(Ref(tree), nodequeue) .+ tree.height

  if any(heights .== Inf)
    error("Inf heights are not supported in plotting")
  end

  tree_x = Vector{Float64}[]
  tree_y = Vector{Float64}[]
  leaf_indices = Int64[]
  node_ids = Int64[]

  for i in eachindex(nodequeue)
    if !isroot(tree, nodequeue[i])
      pn_index = findfirst(parentnode(tree, nodequeue[i]) .== nodequeue)
      push!(tree_x, heights[[i; pn_index; pn_index]])
      push!(tree_y, left_pos[[i; i; pn_index]])
      if isleaf(tree, nodequeue[i])
        push!(leaf_indices, length(tree_y))
      end
      push!(node_ids, nodequeue[i])
    end
  end
  return tree_x, tree_y, leaf_indices, node_ids
 end

@recipe function plot(tree::Tree; plot_leaves = true, hover_ids = false)
  tree_x, tree_y, leaf_indices, node_ids = _treeplot(tree)
  yaxis := false
  grid := :x
  legend := false
  @series begin
    seriestype := :path
    linecolor --> :black
    tree_x, tree_y
  end
  if plot_leaves
    @series begin
      seriestype := :scatter
      marker --> :diamond
      markersize --> min(3, 100 * length(leaf_indices)^-1)
      if hover_ids
        hover := ["Node $x" for x in node_ids[leaf_indices]]
      end
      getindex.(tree_x[leaf_indices], 1),  getindex.(tree_y[leaf_indices], 1)
    end
  end
end


text:
module OctFVM

    using LinearAlgebra
    using SparseArrays
    using IterTools

    include("OctMesh.jl")
    include("AdjacentMats.jl")

end # module


text:
"""User Timeline Endpoint as constant"""
const USER_TIMELINE = "https://api.twitter.com/1.1/statuses/user_timeline.json"

"""
    Authentictor (Datatype:TwiliteTimeline)

Authenticator is defined datatype to create an object which will be used
to authenticate the user when interacting with the twitter API. Someone must
generate the Consumer Key, Consumer Secret, OAuth Token and OAuth Secret
using their own twitter developer account then using those credentionas
create a value of the type Authenticator.

# Example
```
# Here Twitter OAuth token and secret is stored as environment varialbe whith name
# CKEY, CSEC, OTOK, OSEC
julia> creds =  Authentictor(ENV["CKEY"], ENV["CSEC"], ENV["OTOK"], ENV["OSEC"])
```
"""
struct Authentictor
    consumer_key::String
    consumer_secret::String
    oauth_token::String
    oauth_secret::String
end

"""
    ResourceParams(DataType:TwiliteTimeline)

ResourceParams is a defined datatype to store the Resource Parameter mentioned in the twitter API document which we have to send when qurying the API in the request body.

To be able to use twitter API the fild value shoud satisfy the following criteria:
- When Providing valid `count` value the `include_rts` arguments must have to be 1; by default providing any valid `count` value of the `include_rts` will be set to 1 silintly
- When providing constructor argument value for the boolean field such as `trim_user`, `exclude_replies` and `include_rts` the value mush have to be 1 or 0. By the default the values are set to the default value according to Twitter API doc.

# Example
```
julia> resource_params = ResourceParams("realDonaldTrump", count=100)
ResourceParams("realDonaldTrump", nothing, nothing, 100, nothing, false, true)

julia> resource_params = ResourceParams("realDonaldTrump", count=100, exclude_replies=true)
ResourceParams("realDonaldTrump", nothing, nothing, 100, nothing, false, true)
```
"""
struct ResourceParams
    screen_name::String
    user_id::Union{Integer, Nothing}
    since_id::Union{Integer, Nothing}
    count::Union{Integer, Nothing}
    max_id::Union{Integer, Nothing}
    trim_user::Union{Bool, Nothing}
    exclude_replies::Union{Bool, Nothing}

    # Custom Constructor for the Type Signature
    ResourceParams(
        screen_name;
        user_id = nothing,
        science_id = nothing,
        count=1,
        max_id=nothing,
        trim_user=false,
        exclude_replies=true,
    ) = (
    count isa Integer && count < 1
    ? error("Count must have to be greter than or equal to 1")
    : new(screen_name, user_id, science_id, count, max_id, trim_user, exclude_replies)
    )
end

"""
    params_to_dict(rp::ResourceParams)

Is a internal function to convert the user provided Twitter API parameter into dictionary which `oauth_request_resource` from `OAuth` package expect. The method silently add `include_rts` paramater value to 1 when providing an valid `count paramter`.
"""
function params_to_dict(rp::ResourceParams)
    params = Dict{String, Union{String, Integer}}()
    for fn in fieldnames(ResourceParams)
        fv = getfield(rp, fn)
        if fv |> !isnothing
            fn_str = string(fn)
            if fv isa Bool
                params[fn_str] = string(fv)
            else
                params[fn_str] = fv
            end
        end
    end
    if haskey(params, "count")
        params["include_rts"] = 1
    end
    return params
end

"""
    collect_tweets(a::Authentictor, rp::ResourceParams)

Try to collect tweets based on given object of type Autenticator and ResourceParams. The `collect_tweets` is a wrapper around the `oauth_request_resource` method from `OAuth` package which is responsible for sending the `get` request to the `statuses/user_timeline`.

# Example
```
# Here Twitter OAuth token and secret is stored as environment varialbe whith name
# CKEY, CSEC, OTOK, OSEC
julia> creds = Authentictor(ENV["CKEY"], ENV["CSEC"], ENV["OTOK"], ENV["OSEC"])
julia> resource_params = ResourceParams("realDonaldTrump", count=100)
julia> tweets = collect_tweets(creds, resource_params);
```
That should return 100 latest tweets from the profile @realDonaldTrump.
"""
function collect_tweets(a::Authentictor, rp::ResourceParams)
     @warn "use get_tweets method instead, that method will be removed in 1.3 release. Read more details in the readme file."
    # paramter convertion to dictionary
    params = params_to_dict(rp)
    # Response Object
    resp = oauth_request_resource(
        USER_TIMELINE,
        "GET",
        params,
        a.consumer_key,
        a.consumer_secret,
        a.oauth_token,
        a.oauth_secret
    )
    if resp.status == 200
        newdata = JSON.parse(String(resp.body))
        return newdata
    else
        error("Twitter API returned $(r.status) status")
    end
end


text:
using GeodesicRendering
using AccretionGeometry
using CarterBoyerLindquist
using AccretionFormulae
using StaticArrays
using Printf

using Plots
gr()

function temperature(m, gp, max_time; kwargs...)
    g = AccretionFormulae.redshift(m, gp, max_time; kwargs...)

    M = m.M
    a_star = m.a
    temperature = AccretionFormulae.observed_temperature(gp.u[2], a_star, M, g)

end

function radius(m, gp, max_time; kwargs...)
    gp.u[2]
end

function temperature_render(;
    mass = 1,
    spin = 0.998,
    obs_angle = 85.0,
    disc_angle = 90.0,
    tolerance = 1e-8,
    dtmax = 1000.0,
    size_multiplier::Int64 = 1,
    resolution = 400,
    fov = 3.0,
    η = 0.1,
    η_phys = 0.1,
    edd_ratio = 0.1,
    edd_ratio_phys = 0.1,
)

    m = CarterMethodBL(M = 1.0, a = spin)
    M = m.M

    # observer position
    u = [0.0, 1000.0, deg2rad(obs_angle), 0.0]
    R_isco = AccretionFormulae.r_isco(m.a, m.M)

    # disc
    d = GeometricThinDisc(R_isco + 1, 50.0, deg2rad(disc_angle))

    # cache the render
    cache = @time prerendergeodesics(
        m,
        u,
        2000.0,
        d,
        fov_factor = fov * size_multiplier,
        abstol = tolerance,
        reltol = tolerance,
        image_width = 350 * size_multiplier,
        image_height = 250 * size_multiplier,
        dtmax = dtmax,
    )

    # create and compose the PointFunctions
    temperature_vf = (
        PointFunction(temperature) ∘
        FilterPointFunction((m, gp, max_time; kwargs...) -> gp.u[2] > R_isco, NaN) ∘
        ConstPointFunctions.filter_early_term
    )

    radius_vf = (
        PointFunction(radius) ∘
        FilterPointFunction((m, gp, max_time; kwargs...) -> gp.u[2] > R_isco, NaN) ∘
        ConstPointFunctions.filter_early_term
    )

    radius_img = GeodesicRendering.apply(radius_vf, cache)
    temperature_img = GeodesicRendering.apply(temperature_vf, cache)

    # correcting for physical mass
    M_phys = mass * 1.99e30
    r_isco_phys = AccretionFormulae.r_isco(M_phys, 0.998)
    r_g_phys = M_phys

    numerators =
        AccretionFormulae.mass_scale_fraction.(
            M_phys,
            η_phys,
            edd_ratio_phys,
            r_isco_phys,
            r_g_phys,
            radius_img * M_phys,
        )
    denominators =
        AccretionFormulae.mass_scale_fraction.(M, η, edd_ratio, R_isco, M, radius_img)

    fractions = numerators ./ denominators
    temperature_img .*= fractions

    # # autoscale
    # scale = maximum(filter(!isnan,temperature_img))
    # scale = floor(log(10, scale))
    # scale = 10^scale

    # scaling image
    scale = 1e7
    scalestr = @sprintf "%.E" scale
    new_img = reverse(temperature_img, dims = 1)
    new_img ./= scale

    heatmap(
        new_img,
        aspect_ratio = 1.0,
        size = (resolution * 3 / 2, resolution),
        clim = (0, 3),
    )
    # contour(new_img, aspect_ratio=1.0, size=(resolution*3/2, resolution), clim=(0,3))
    title!("Temperature Scale = $scalestr, Mass = $mass M_☼")
end

temperature_render(obs_angle = 85.0, mass = 1)


text:
module AbstractLattices

export ∧, ∨, dist, wedge, vee

function wedge end
function vee end

const ∧ = wedge
const ∨ = vee

function dist end

end


text:
# Central coefficient
ccentral(coeff :: VoF) = 1 / (length(coeff) + 1)

# Sum of all filter coefficients
csum(coeff :: VoF) = ccentral(coeff) + 2sum(coeff)

# Filter evaluated at point `z`
fn(coeff :: VoF, z :: Complex) =
    ccentral(coeff) + sum(c*(z^n + z^-n) for (n, c) in enumerate(coeff))

cnonnegative(coeff :: VoF) =
    2000 * max(0, -minimum(real(fn(coeff, exp(2π * x * im))) for x in 0:0.001:1))

# The target function for coefficients of our filter incorporates the
# following:
#
# 1) We want the sum of all coefficients to be 1 so our filter is
#    low-pass.
# 2) The filter also has to satisfy conditions of Riesz-Fejér
#    factorisation lemma, i.e. map unit circle to non-negative real
#    values.
# 3) The free coefficient is fixed at 1/length because porosity is
#    everything.
#
# The coefficients found by minimizing `targetfn` are not unique.
function targetfn(coeff :: VoF)
    return cnonnegative(coeff) + (csum(coeff) - 1)^2
end

"""
    filter_coeffs(length[; initial])

Get coefficients of symmetrical low-pass FIR filter of length
`2length-1` which maps autocorrelation functions to autocorrelation
functions. Argument `length` being equal to `n` is enough for
upscaling autocorrelations function `n` times.

A keyword argument `initial` can be an array of length `length - 1`
which works as an initial guess for filter coefficients.
"""
function filter_coeffs(length  :: Integer;
                       initial :: VoF = fill(0.95/(2length-1), length-1))
    res = optimize(targetfn, initial)
    return Optim.minimizer(res)
end

# Filter length function
filter_length(ratio) = max(numerator(ratio), denominator(ratio))


text:
using SparseArrays
using ViewerGL
GL = ViewerGL
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

# input generation
V,EV = Lar.randomcuboids(10, .35)
V = GL.normalize2(V,flag=true)
VV = [[k] for k=1:size(V,2)]
GL.VIEW( GL.numbering(.05)((V,[VV, EV]),GL.COLORS[1]) )
@show V
@show EV
# subdivision of input edges
W = convert(Lar.Points, V')
cop_EV = Lar.coboundary_0(EV::Lar.Cells)
cop_EW = convert(Lar.ChainOp, cop_EV)
V, copEV, copFE = Lar.planar_arrangement(W::Lar.Points, cop_EW::Lar.ChainOp)

# compute containment graph of components
bicon_comps = Lar.Arrangement.biconnected_components(copEV)

# visualization of component graphs
EW = Lar.cop2lar(copEV)
W = convert(Lar.Points, V')
comps = [ GL.GLLines(W,EW[comp],GL.COLORS[(k-1)%12+1]) for (k,comp) in enumerate(bicon_comps) ]
GL.VIEW(comps);

# visualization of numbered arrangement
VV = [[k] for k=1:size(V,1)]
EV = Lar.cop2lar(copEV);
FE = Lar.cop2lar(copFE);
FV = cat([[EV[e] for e in face] for face in FE])

# final solid visualization
triangulated_faces = Lar.triangulate2D(V, [copEV, copFE])
V = convert(Lar.Points, V')
FVs = convert(Array{Lar.Cells}, triangulated_faces)
GL.VIEW(GL.GLExplode(V,FVs,1.2,1.2,1.2,99,1));

# polygonal face boundaries
EVs = Lar.FV2EVs(copEV, copFE)
EVs = convert(Array{Array{Array{Int64,1},1},1}, EVs)
GL.VIEW(GL.GLExplode(V,EVs,1.2,1.2,1.2,1,1));


text:
n = 2039;
n1 = n;
# wav = wavs[2]; k = ks[1]
@testset "Delta Spikes" begin
    wavs = (ContinuousWavelets.wavelet(cDb2), ContinuousWavelets.wavelet(cCoif4), ContinuousWavelets.wavelet(cBeyl))
    ks = (1093, 408)
    @testset "at $k, with type $(wav.waveType)" for k in ks, wav in wavs
        x = zeros(n)
        x[k] = 1
        res = 3
        with_logger(ConsoleLogger(stderr, Logging.Error)) do
            res = ContinuousWavelets.cwt(x, wav)
        end
        peaks = argmax(abs.(res), dims = 1)
        @test k - 2 <= peaks[end][1] <= k + 2 # give slight range to handle rounding issues
        Ŵ, ω = (3, 1)
        with_logger(ConsoleLogger(stderr, Logging.Error)) do
            Ŵ, ω = ContinuousWavelets.computeWavelets(n, wav)
        end
        W = irfft(Ŵ, 2n, 1)
        @test !any(abs.(circshift(W[:, 3], k - 1)[1:n] .- res[:, 3]) .> 1e-10)
        # make sure the wavelets don't have support out into the mirrored signal
    end
end


text:
using Test
using YaoAPI
using MLStyle
using YaoLocations
using YaoCompiler
using GPUCompiler
using YaoCompiler.Intrinsics
using CompilerPluginTools
using YaoCompiler.Intrinsics: measure

@device function test_basic(theta, phi)
    # syntax sugar
    1 => X
    @gate 2 => Z
    @ctrl 1 4 => Rx(theta)
    @ctrl 2 4 => Ry(phi)
    a = @measure 3
    # direct intrinsic
    apply(Y, 3)
    apply(X, 1, 4)
    c = measure(2)
    return (a = a, b = c)
end

@device function test_pure_quantum()
    ret = @gate 1:4 => test_basic(1.0, 2.0)
    @ctrl 2 1 => Rx(2.2)
    return ret
end

ci, type = @yao_code_typed(test_pure_quantum())[1]

@device function routine2(theta)
    2 => Rx(theta)
end

ci, type = @yao_code_typed(routine2(2.0))[1]

op = test_pure_quantum()
interp = YaoInterpreter()
ci, type = @yao_code_typed(op)[1]

@test_codeinfo ci begin
    Expr(
        :invoke,
        _,
        GlobalRef(Intrinsics, :apply),
        Argument(2),
        QuoteNode(X),
        QuoteNode(Locations(1)),
    )::Nothing
    Expr(
        :invoke,
        _,
        GlobalRef(Intrinsics, :apply),
        Argument(2),
        QuoteNode(Z),
        QuoteNode(Locations(2)),
    )::Nothing
    Expr(
        :invoke,
        _,
        GlobalRef(Intrinsics, :apply),
        Argument(2),
        QuoteNode(Rx(1.0)),
        QuoteNode(Locations(4)),
        QuoteNode(&(CtrlLocations(1))),
    )::Nothing
    Expr(
        :invoke,
        _,
        GlobalRef(Intrinsics, :apply),
        Argument(2),
        QuoteNode(Ry(2.0)),
        QuoteNode(Locations(4)),
        QuoteNode(&(CtrlLocations(2))),
    )::Nothing
    Expr(
        :invoke,
        _,
        GlobalRef(Intrinsics, :measure),
        Argument(2),
        QuoteNode(Locations(3)),
    )::MeasureResult{Int}
    Expr(
        :invoke,
        _,
        GlobalRef(Intrinsics, :apply),
        Argument(2),
        QuoteNode(Y),
        QuoteNode(Locations(3)),
    )::Nothing
end

struct JLDummyTarget <: YaoCompileTarget end

function test_pass()
    @test true
    return
end

function YaoCompiler.target_specific_pipeline(::JLDummyTarget, ir::IRCode)
    mi = method_instances(test_pass, ())[1]
    for i in 1:length(ir.stmts)
        e = ir.stmts[i][:inst]
        @switch e begin
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :measure), reg, locs)
            ir.stmts[i][:inst] = QuoteNode(MeasureResult(5))
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :barrier), reg, locs)
            ir.stmts[i][:inst] = Expr(:invoke, mi, GlobalRef(Main, :test_pass))
            ir.stmts[i][:type] = Nothing
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), reg, gate, locs)
            ir.stmts[i][:inst] = Expr(:invoke, mi, GlobalRef(Main, :test_pass))
            ir.stmts[i][:type] = Nothing
            @case Expr(:invoke, _, GlobalRef(Intrinsics, :apply), reg, gate, locs, ctrl)
            ir.stmts[i][:inst] = Expr(:invoke, mi, GlobalRef(Main, :test_pass))
            ir.stmts[i][:type] = Nothing
            @case _
            nothing
        end
    end
    return ir
end

f = YaoCompiler.compile(JLDummyTarget(), Intrinsics.apply, Tuple{AnyReg,typeof(op)})
test_ret = f(AnyReg(), op)
@test test_ret.a.result == 5
@test test_ret.b.result == 5


text:
module InvHessian

using Test
using DiffPrecTest
using ProximalBase
using CovSel
using CoordinateDescent
using Random, LinearAlgebra, SparseArrays, Statistics, Distributions

function genData(p)
  Sigmax = Matrix(1.0I, p, p)
  Sigmay = zeros(p,p)
  rho = 0.7
  for i=1:p
    for j=1:p
      Sigmay[i,j]=rho^abs(i-j)
    end
  end
  sqmy = sqrt(Sigmay)
  n = 1000
  X = randn(n,p)
  Y = randn(n,p) * sqmy;
  hSx = cov(X)
  hSy = cov(Y)
  hSx, hSy
end

function Q(p::Int)
  _c = 1. / sqrt(2.)
  rp = div((p+1)*p,2)
  Q = spzeros( rp,  p*p )
  ind = 0
  for ci=1:p
    for ri=ci:p
      ind += 1
      if ri == ci
        Q[ind, (ci-1)*p + ri] = 1.
      else
        Q[ind, (ci-1)*p + ri] = _c
        Q[ind, (ri-1)*p + ci] = _c
      end
    end
  end
  Q
end

function Qt(p::Int)
  _c = 1. / sqrt(2.)
  rp = div((p+1)*p,2)
  Q = spzeros( p*p, rp )
  ind = 0
  for ci=1:p
    for ri=ci:p
      ind += 1
      if ri == ci
        Q[(ci-1)*p + ri, ind] = 1.
      else
        Q[(ci-1)*p + ri, ind] = _c
        Q[(ri-1)*p + ci, ind] = _c
      end
    end
  end
  Q
end

@testset "invert Kroneker" begin
  for rep=1:50
    p = 20
    hSx, hSy = genData(p)
    A = kron(hSy, hSx)

    λ = rand(Uniform(0.05, 0.3))
    g = ProximalBase.ProxL1(λ)

    for i = [(1,1), (3,1), (4,5)]
      ri = i[1]
      ci = i[2]
      x = ProximalBase.SparseIterate(p*p)
      x1 = ProximalBase.SparseIterate(p*p)

      f = DiffPrecTest.CDInverseKroneckerLoss(Symmetric(hSx), Symmetric(hSy), ri, ci)
      b = zeros(Float64, p*p)
      b[(ci-1)*p+ri] = -1.

      f1 = CoordinateDescent.CDQuadraticLoss(A, b)
      CoordinateDescent.coordinateDescent!(x, f, g, CoordinateDescent.CDOptions(;maxIter=5000, optTol=1e-12))
      CoordinateDescent.coordinateDescent!(x1, f1, g, CoordinateDescent.CDOptions(;maxIter=5000, optTol=1e-12))

      @test convert(Vector, x) ≈ convert(Vector, x1) atol=1e-7
    end
  end
end
#
# @testset "invert SymKroneker" begin
#   for rep=1:50
#     p = 20
#     hSx, hSy = genData(p)
#     A = (kron(hSy, hSx) + kron(hSx, hSy)) / 2.
#
#     λ = rand(Uniform(0.05, 0.3))
#     g = ProximalBase.ProxL1(λ)
#
#     for i = [(1,1), (3,1), (4,5)]
#       ri = i[1]
#       ci = i[2]
#       x = ProximalBase.SparseIterate(p*p)
#       x1 = ProximalBase.SparseIterate(p*p)
#
#       f = DiffPrecTest.CDInverseSymKroneckerLoss(Symmetric(hSx), Symmetric(hSy), ri, ci)
#       b = zeros(Float64, p*p)
#       b[(ci-1)*p+ri] = -1.
#
#       f1 = CoordinateDescent.CDQuadraticLoss(A, b)
#       CoordinateDescent.coordinateDescent!(x, f, g, CoordinateDescent.CDOptions(;maxIter=5000, optTol=1e-12))
#       CoordinateDescent.coordinateDescent!(x1, f1, g, CoordinateDescent.CDOptions(;maxIter=5000, optTol=1e-12))
#
#       @test convert(Vector, x) ≈ convert(Vector, x1) atol=1e-7
#     end
#   end
# end
#
# @testset "invert reduced SymKroneker" begin
#   for rep=1:50
#     p = 20
#     hSx, hSy = genData(p)
#     A = (kron(hSy, hSx) + kron(hSx, hSy)) / 2.
#
#     λ = rand(Uniform(0.05, 0.3))
#     g = ProximalBase.ProxL1(λ)
#
#     for i = [(1,1), (3,1), (4,5)]
#       ri = i[1]
#       ci = i[2]
#       x = ProximalBase.SparseIterate(p*p)
#       x1 = ProximalBase.SparseIterate(p*p)
#
#       f = DiffPrecTest.CDInverseSymKroneckerLoss(Symmetric(hSx), Symmetric(hSy), ri, ci)
#       b = zeros(Float64, p*p)
#       b[(ci-1)*p+ri] = -1.
#
#       f1 = CoordinateDescent.CDQuadraticLoss(A, b)
#       CoordinateDescent.coordinateDescent!(x, f, g, CoordinateDescent.CDOptions(;maxIter=5000, optTol=1e-12))
#       CoordinateDescent.coordinateDescent!(x1, f1, g, CoordinateDescent.CDOptions(;maxIter=5000, optTol=1e-12))
#
#       @test convert(Vector, x) ≈ convert(Vector, x1) atol=1e-7
#     end
#   end
# end

end


text:
abstract GP
include("full_gp.jl")
include("laplacian_gp.jl")
include("nystrom_gp.jl")
include("fitc_gp.jl")


text:
using XRayQuadrics
using Plots
using Colors

# macro expands to:
#   const PlotTQ{ArgTypes} = Combined{plottq, ArgTypes}
#   plottq(args...; kw_args...) = ...
#   plottq!(args...; kw_args...) = ...
#   function default_theme(scene, ::PlotTQ) = ...

# @Makie.recipe(PlotTQ, obj) do scene
#     Makie.Attributes(;
#         plot_caps = true,
#         color = :red,
#         seriesalpha = 0.667,
#         light_source = [0;0;1]
#     )
# end

# Makie.plottype(::TruncatedQuadric) = PlotTQ{<:Tuple{TruncatedQuadric}}

# function Makie.plot!(plot::PlotTQ{<:Tuple{TruncatedQuadric}})
#     tq = plot[:obj][]

#     (Xs, Ys, Zs) = cartesian_grid(tq)

#     X = Xs[1]
#     Y = Ys[1]
#     Z = Zs[1]

#     Makie.surface!(
#         plot, X, Y, Z,
#         shading=false,
#         color=fill(RGBA(1.0,0.0,0.0,0.4), size(X)...)
#     )
#     if plot[:plot_caps][]
#         Makie.surface!(
#             plot, Xs[2], Ys[2], Zs[2],
#             shading=false,
#             color=fill(RGBA(0.0,0.0,1.0,0.4), size(Xs[2]))
#         )
#         Makie.surface!(
#             plot, Xs[3], Ys[3], Zs[3],
#             shading=false,
#             color=fill(RGBA(0.0,0.0,1.0,0.4), size(Xs[3]))
#         )
#     end
#     return plot
# end

Plots.@recipe function f(tq::TruncatedQuadric; plot_caps = false)
    # color       --> :blue
    seriesalpha --> 0.667
    seriestype  :=  :surface
    legend      --> false
    # markershape --> (add_marker ? :circle : :none)
    delete!(plotattributes, :add_marker)
    (X, Y, Z) = cartesian_grid(tq)
    @series begin
        seriestype  := :surface
        legend      --> false
        X[1], Y[1], Z[1]
    end
    # if plot_caps
        @series begin
            seriestype  := :surface
            legend      --> false
            X[2], Y[2], Z[2]
        end
        @series begin
            seriestype  := :surface
            legend      --> false
            X[3], Y[3], Z[3]
        end
    # end
end

function cartesian_grid(tq::TruncatedQuadric)
    q = tq.q
    ps = tq.p
    s = changerepresentation(q)
    T = typeof(s)

    (nθ, nζ) = (30, 100)
    (Xs, Ys, Zs) = get_mesh(s, ps, nθ, nζ)
    X = Xs[1]
    Y = Ys[1]
    Z = Zs[1]
    return (Xs, Ys, Zs)
end

function get_mesh(s::Union{Cylinder, Paraboloid, Hyperboloid}, ps::Vector{Plane}, nθ::Int64, nζ::Int64)

    ca1 = axis_plane_intersection(s, ps[1])
    ca2 = axis_plane_intersection(s, ps[2])

    h = norm(ca2 - ca1)

    h1 = ca1'*s.a
    h2 = ca2'*s.a

    a = s.a
    c = s.c
    θ = range(0, stop=2π, length=nθ)
    # ζ = range(0, stop=h, length=nζ)
    ζ = range(h1, stop=h2, length=nζ)

    X = zeros(nθ, nζ)
    Y = zeros(nθ, nζ)
    Z = zeros(nθ, nζ)

    if typeof(s) == Cylinder
        R = s.R
        for i in 1:nθ
            for j in 1:nζ
                X[i,j] = R*cos(θ[i])
                Y[i,j] = R*sin(θ[i])
                Z[i,j] = ζ[j]
            end
        end
    elseif typeof(s) == Paraboloid
        b = s.b
        for i in 1:nθ
            for j in 1:nζ
                X[i,j] = b*sqrt(ζ[j])*cos(θ[i])
                Y[i,j] = b*sqrt(ζ[j])*sin(θ[i])
                Z[i,j] = ζ[j]
            end
        end
    elseif typeof(s) == Hyperboloid
        R = s.R
        b = s.b
        for i in 1:nθ
            for j in 1:nζ
                X[i,j] = R*sqrt(1 + ζ[j]^2/b^2)*cos(θ[i])
                Y[i,j] = R*sqrt(1 + ζ[j]^2/b^2)*sin(θ[i])
                Z[i,j] = ζ[j]
            end
        end
    end

    Xc1 = [zeros(nθ, 1) X[:,1]]
    Xc2 = [zeros(nθ, 1) X[:,end]]
    Yc1 = [zeros(nθ, 1) Y[:,1]]
    Yc2 = [zeros(nθ, 1) Y[:,end]]
    Zc1 = [Z[:,1] Z[:,1]]
    Zc2 = [Z[:,1] Z[:,1]]

    # print(typeof(s))
    (X, Y, Z) = transform_to_axis(X, Y, Z, [0;0;1], a, c)
    (Xc1, Yc1, Zc1) = transform_to_axis(Xc1, Yc1, Zc1, [0;0;1], ps[1].a, c)   # FIX THESE: WRONG MESH POSITION
    (Xc2, Yc2, Zc2) = transform_to_axis(Xc2, Yc2, Zc2, [0;0;1], ps[2].a, c + h2*a)

    return ((X,Xc1,Xc2), (Y,Yc1,Yc2), (Z,Zc1,Zc2))
end

function axis_plane_intersection(s::Union{Cylinder, Paraboloid, Hyperboloid}, p::Plane)
    # project plane points onto axis
    ca = s.c + (p.a'*p.c - p.a'*s.c)/(p.a'*s.a)*s.a
    return ca
end

function transform_to_axis(X, Y, Z, oldax, newax, center)
    # can do checks on sizes of X, Y, Z
    (nθ, nζ) = size(X)
    transform = zeros(3,3)
    if newax'*oldax == 1
        transform = I
    elseif newax'*oldax == -1
        transform = [1  0  0;
                     0  0 -1;
                     0 -1  0]
    else
        # cylinder axis oblique to z-axis:
        ax = cross(oldax, newax)
        ax = ax./norm(ax)
        cosang = newax'*oldax
        # Rodrigues's formula/axis-angle representation:
        axso3 = [ 0     -ax[3]   ax[2];
                  ax[3]  0      -ax[1];
                 -ax[2]  ax[1]   0]
        transform = I*cosang + (1 - cosang)*(ax*ax') + axso3*sqrt(1 - cosang^2)
    end

    # transform all cylinder points:
    for i = 1:nθ
        for j = 1:nζ
            newpos = transform*[X[i,j]; Y[i,j]; Z[i,j]]
            X[i,j] = newpos[1] + center[1]
            Y[i,j] = newpos[2] + center[2]
            Z[i,j] = newpos[3] + center[3]
        end
    end
    return (X, Y, Z)
end


text:
"""
All steps of data cleaning that most researchers do can be performed using the conventional cleaning funciton.

# Example
```julia
radiance_datacube = rand(1:1000, 10, 10, 10)
clouds_datacube = rand(1:1000, 10, 10, 10)
conventional_cleaning(radiance_datacube, clouds_datacube)
```
"""
function conventional_cleaning(radiance_datacube, clouds_datacube)
    tmp = mark_missing(radiance_datacube, clouds_datacube)
    GC.gc()
    tmp =  replace_negative(tmp)
    GC.gc()
    noise = background_noise_mask(tmp, clouds_datacube, 0.4)
    GC.gc()
    tmp = apply_mask(tmp, noise)
    GC.gc()
    stable_pixels = outlier_mask(tmp, noise)
    GC.gc()
    tmp = apply_mask(tmp, stable_pixels)
    GC.gc()
    tmp = long_apply(outlier_ts, tmp)
    GC.gc()
    tmp = long_apply(linear_interpolation, tmp)
    GC.gc()
    return Array{Float16}(tmp)
end
"""
The PatnaikSTT2021 function performs all the steps of the new cleaning procedure described in [But clouds got in my way: Bias and bias correction of VIIRS nighttime lights data in the presence of clouds, Ayush Patnaik, Ajay Shah, Anshul Tayal, Susan Thomas](https://www.xkdr.org/releases/PatnaikShahTayalThomas_2021_bias_correction_nighttime_lights.html) as conventional cleaning.

# Example
```julia
radiance_datacube = rand(1:1000, 10, 10, 10)
clouds_datacube = rand(1:1000, 10, 10, 10)
PatnaikSTT2021(radiance_datacube, clouds_datacube)
```
"""
function PatnaikSTT2021(radiance_datacube, clouds_datacube)
    tmp = mark_missing(radiance_datacube, clouds_datacube)
    GC.gc()
    tmp = replace_negative(tmp)
    GC.gc()
    noise = background_noise_mask(tmp, clouds_datacube, 0.4)
    GC.gc()
    tmp = apply_mask(tmp, noise)
    GC.gc()
    stable_pixels = outlier_mask(tmp, noise)
    GC.gc()
    tmp = apply_mask(tmp, stable_pixels)
    GC.gc()
    tmp = long_apply(outlier_ts, tmp)
    GC.gc()
    mask = noise .* stable_pixels
    GC.gc()
    tmp = bias_correction(tmp, clouds_datacube, mask)
    GC.gc()
    tmp = long_apply(linear_interpolation, tmp)
    GC.gc()
    return Array{Float16}(tmp)
end


text:
using Test
using PETScX

@testset "init tests" begin
    for petsclib in PETScX.petsclibs
        # The first time through finalize should be false since we have never
        # initialized petsclib yet...
        initial_finalized_value = false

        for x in (petsclib, typeof(petsclib))
            # since we haven't called anything these should be false!
            @test !(PETScX.Initialized(x))
            @test PETScX.Finalized(x) == initial_finalized_value

            # The second time through  time through finalize should be true since
            # we have initialized petsclib yet...
            initial_finalized_value = true

            # Initialize PETSc
            PETScX.Initialize(x)

            # Check values again
            @test PETScX.Initialized(x)
            @test !(PETScX.Finalized(x))

            PETScX.Finalize(x)

            # Check values again
            @test !(PETScX.Initialized(x))
            @test PETScX.Finalized(x)
        end
    end
end


text:
## types for  diffusionmap calculation

struct Diffusionmap
    data::Matrix
    kernel::AbstractKernel
    laplace_type::AbstractLaplacian
    threshold::Int64
end

Diffusionmap(data; kernel::AbstractKernel=InverseDistanceKernel(), laplace_type::AbstractLaplacian=RowNormalizedLaplacian(), threshold::Int64=size(data,1)) = Diffusionmap(data, kernel, laplace_type, threshold)
Diffusionmap(data, func::Function; laplace_type::AbstractLaplacian=RowNormalizedLaplacian(), threshold::Int64=size(data,1)) = Diffusionmap(data, CustomKernel(func), laplace_type, threshold)

function Base.show(io::IO, ::MIME"text/plain", dp::Diffusionmap)
    println(io, "Diffusionmap")
    println(io, "Kernel = $(typeof(dp.kernel))")
    println(io, "Laplace type = $(dp.laplace_type)")
    println(io, "Threshold = $(dp.threshold)")
end


text:
int_rules_1_1_4_2 = @theory begin

    #= ::Subsection::Closed:: =#
    #= 1.1.4.2*(c*x)^m*(a*x^j+b*x^n)^p =#
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (1 / ~n) * Subst(Antiderivative((~a * (~x) ^ Simplify(~j / ~n) + ~b * ~x) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && EqQ(Simplify((~m - ~n) + 1), 0))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => (-((~c) ^ (~j - 1)) * (~c * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~a * (~n - ~j) * (~p + 1)) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (EqQ(((~m + ~n * ~p + ~n) - ~j) + 1, 0) && (IntegerQ(~j) || GtQ(~c, 0)))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => (-((~c) ^ (~j - 1)) * (~c * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~a * (~n - ~j) * (~p + 1)) + (((~c) ^ ~j * (((~m + ~n * ~p + ~n) - ~j) + 1)) / (~a * (~n - ~j) * (~p + 1))) * Antiderivative((~c * ~x) ^ (~m - ~j) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (ILtQ(Simplify((((~m + ~n * ~p + ~n) - ~j) + 1) / (~n - ~j)), 0) && (LtQ(~p, -1) && (IntegerQ(~j) || GtQ(~c, 0))))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => ((~c) ^ (~j - 1) * (~c * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~a * (~m + ~j * ~p + 1)) - ((~b * (((~m + ~n * ~p + ~n) - ~j) + 1)) / (~a * (~c) ^ (~n - ~j) * (~m + ~j * ~p + 1))) * Antiderivative((~c * ~x) ^ ((~m + ~n) - ~j) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (ILtQ(Simplify((((~m + ~n * ~p + ~n) - ~j) + 1) / (~n - ~j)), 0) && (NeQ(~m + ~j * ~p + 1, 0) && (IntegersQ(~j, ~n) || GtQ(~c, 0))))))
    @apply_utils Antiderivative((~c * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => (((~c) ^ IntPart(~m) * (~c * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && ILtQ(Simplify((((~m + ~n * ~p + ~n) - ~j) + 1) / (~n - ~j)), 0)))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (1 / ~n) * Subst(Antiderivative((~x) ^ (Simplify((~m + 1) / ~n) - 1) * (~a * (~x) ^ Simplify(~j / ~n) + ~b * ~x) ^ ~p, ~x), ~x, (~x) ^ ~n) <-- FreeQ([~a, ~b, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify((~m + 1) / ~n)) && NeQ((~n) ^ 2, 1)))))
    @apply_utils Antiderivative((~c * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~c) ^ IntPart(~m) * (~c * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && (IntegerQ(Simplify((~m + 1) / ~n)) && NeQ((~n) ^ 2, 1)))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~m * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => ((~c * ~x) ^ (~m + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p) / (~c * (~m + ~j * ~p + 1)) - ((~b * ~p * (~n - ~j)) / ((~c) ^ ~n * (~m + ~j * ~p + 1))) * Antiderivative((~c * ~x) ^ (~m + ~n) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (Not(IntegerQ(~p)) && (LtQ(0, ~j, ~n) && ((IntegersQ(~j, ~n) || GtQ(~c, 0)) && (GtQ(~p, 0) && LtQ(~m + ~j * ~p + 1, 0)))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => ((~c * ~x) ^ (~m + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p) / (~c * (~m + ~n * ~p + 1)) + ((~a * (~n - ~j) * ~p) / ((~c) ^ ~j * (~m + ~n * ~p + 1))) * Antiderivative((~c * ~x) ^ (~m + ~j) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (Not(IntegerQ(~p)) && (LtQ(0, ~j, ~n) && ((IntegersQ(~j, ~n) || GtQ(~c, 0)) && (GtQ(~p, 0) && NeQ(~m + ~n * ~p + 1, 0)))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => ((~c) ^ (~n - 1) * (~c * ~x) ^ ((~m - ~n) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * (~n - ~j) * (~p + 1)) - (((~c) ^ ~n * (((~m + ~j * ~p) - ~n) + ~j + 1)) / (~b * (~n - ~j) * (~p + 1))) * Antiderivative((~c * ~x) ^ (~m - ~n) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c], ~x) && (Not(IntegerQ(~p)) && (LtQ(0, ~j, ~n) && ((IntegersQ(~j, ~n) || GtQ(~c, 0)) && (LtQ(~p, -1) && GtQ(~m + ~j * ~p + 1, ~n - ~j)))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => (-((~c) ^ (~j - 1)) * (~c * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~a * (~n - ~j) * (~p + 1)) + (((~c) ^ ~j * (((~m + ~n * ~p + ~n) - ~j) + 1)) / (~a * (~n - ~j) * (~p + 1))) * Antiderivative((~c * ~x) ^ (~m - ~j) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~m], ~x) && (Not(IntegerQ(~p)) && (LtQ(0, ~j, ~n) && ((IntegersQ(~j, ~n) || GtQ(~c, 0)) && LtQ(~p, -1))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => ((~c) ^ (~n - 1) * (~c * ~x) ^ ((~m - ~n) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~b * (~m + ~n * ~p + 1)) - ((~a * (~c) ^ (~n - ~j) * (((~m + ~j * ~p) - ~n) + ~j + 1)) / (~b * (~m + ~n * ~p + 1))) * Antiderivative((~c * ~x) ^ (~m - (~n - ~j)) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (Not(IntegerQ(~p)) && (LtQ(0, ~j, ~n) && ((IntegersQ(~j, ~n) || GtQ(~c, 0)) && (GtQ(((~m + ~j * ~p + 1) - ~n) + ~j, 0) && NeQ(~m + ~n * ~p + 1, 0)))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => ((~c) ^ (~j - 1) * (~c * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~a * (~m + ~j * ~p + 1)) - ((~b * (((~m + ~n * ~p + ~n) - ~j) + 1)) / (~a * (~c) ^ (~n - ~j) * (~m + ~j * ~p + 1))) * Antiderivative((~c * ~x) ^ ((~m + ~n) - ~j) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~m, ~p], ~x) && (Not(IntegerQ(~p)) && (LtQ(0, ~j, ~n) && ((IntegersQ(~j, ~n) || GtQ(~c, 0)) && LtQ(~m + ~j * ~p + 1, 0))))
    @apply_utils Antiderivative((~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (1 / (~m + 1)) * Subst(Antiderivative((~a * (~x) ^ Simplify(~j / (~m + 1)) + ~b * (~x) ^ Simplify(~n / (~m + 1))) ^ ~p, ~x), ~x, (~x) ^ (~m + 1)) <-- FreeQ([~a, ~b, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && (NeQ(~m, -1) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n)))))))
    @apply_utils Antiderivative((~c * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~c) ^ IntPart(~m) * (~c * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && (IntegerQ(Simplify(~j / ~n)) && (NeQ(~m, -1) && (IntegerQ(Simplify(~n / (~m + 1))) && Not(IntegerQ(~n)))))))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => ((~c * ~x) ^ (~m + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p) / (~c * ~p * (~n - ~j)) + (~a / (~c) ^ ~j) * Antiderivative((~c * ~x) ^ (~m + ~j) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p - 1), ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n], ~x) && (IGtQ(~p + 1 / 2, 0) && (NeQ(~n, ~j) && (EqQ(Simplify(~m + ~j * ~p + 1), 0) && (IntegerQ(~j) || GtQ(~c, 0)))))
    @apply_utils Antiderivative((~x) ^ ~(m') / sqrt(~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')), ~x) => (-2 / (~n - ~j)) * Subst(Antiderivative(1 / (1 - ~a * (~x) ^ 2), ~x), ~x, (~x) ^ (~j / 2) / sqrt(~a * (~x) ^ ~j + ~b * (~x) ^ ~n)) <-- FreeQ([~a, ~b, ~j, ~n], ~x) && (EqQ(~m, ~j / 2 - 1) && NeQ(~n, ~j))
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => (-((~c) ^ (~j - 1)) * (~c * ~x) ^ ((~m - ~j) + 1) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1)) / (~a * (~n - ~j) * (~p + 1)) + (((~c) ^ ~j * (((~m + ~n * ~p + ~n) - ~j) + 1)) / (~a * (~n - ~j) * (~p + 1))) * Antiderivative((~c * ~x) ^ (~m - ~j) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ (~p + 1), ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n], ~x) && (ILtQ(~p + 1 / 2, 0) && (NeQ(~n, ~j) && (EqQ(Simplify(~m + ~j * ~p + 1), 0) && (IntegerQ(~j) || GtQ(~c, 0)))))
    @apply_utils Antiderivative((~c * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~n) ^ ~p, ~x) => (((~c) ^ IntPart(~m) * (~c * ~x) ^ FracPart(~m)) / (~x) ^ FracPart(~m)) * Antiderivative((~x) ^ ~m * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (IntegerQ(~p + 1 / 2) && (NeQ(~n, ~j) && EqQ(Simplify(~m + ~j * ~p + 1), 0)))
    #= Antiderivative((~x)^(~m')*((~a')*(~x)^(~j')+(~b')*(~x)^(~n'))^(~p),~x) :=  (a*x^j+b*x^n)^(p+1)/(b*p*(n-j)*x^(n+j*p))*HypergeometricFunctions._₂F₁(1,1,1-p,- a/(b*x^(n-j))) <-- FreeQ([a,b,j,m,n,p],x) && NeQ(n,j) &&  EqQ(m+j*p+1,0) =#
    #= Antiderivative((~x)^(~m')*((~a')*(~x)^(~j')+(~b')*(~x)^(~n'))^(~p),~x) :=  (a*x^j+b*x^n)^(p+1)/(b*(p-1)*(n-j)*x^(2*n+j*(p-1)))*HypergeometricFunctions._₂F₁( 1,2,2-p,-a/(b*x^(n-j))) <-- FreeQ([a,b,j,m,n,p],x) && NeQ(n,j) &&  EqQ(m+n+(p-1)*j+1,0) =#
    #= Antiderivative((~x)^(~m')*((~a')*(~x)^(~j')+(~b')*(~x)^(~n'))^(~p),~x) :=  (x^(m-j+1)*(a*x^j+b*x^n)^(p+1))/(a*(m+j*p+1))*HypergeometricFunctions._₂F₁(1,(m+ n*p+1)/(n-j)+1,(m+j*p+1)/(n-j)+1,-b*x^(n-j)/a) <--  FreeQ([a,b,j,m,n,p],x) && NeQ(n,j) && NeQ(m+j*p+1,0) &&  NeQ(m+n+(p-1)*j+1,0) =#
    @apply_utils Antiderivative((~(c') * ~x) ^ ~(m') * (~(a') * (~x) ^ ~(j') + ~(b') * (~x) ^ ~(n')) ^ ~p, ~x) => (((~c) ^ IntPart(~m) * (~c * ~x) ^ FracPart(~m) * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ FracPart(~p)) / ((~x) ^ (FracPart(~m) + ~j * FracPart(~p)) * (~a + ~b * (~x) ^ (~n - ~j)) ^ FracPart(~p))) * Antiderivative((~x) ^ (~m + ~j * ~p) * (~a + ~b * (~x) ^ (~n - ~j)) ^ ~p, ~x) <-- FreeQ([~a, ~b, ~c, ~j, ~m, ~n, ~p], ~x) && (Not(IntegerQ(~p)) && (NeQ(~n, ~j) && PosQ(~n - ~j)))
    @apply_utils Antiderivative((~u) ^ ~(m') * (~(a') * (~v) ^ ~(j') + ~(b') * (~v) ^ ~(n')) ^ ~(p'), ~x) => ((~u) ^ ~m / (Coefficient(~v, ~x, 1) * (~v) ^ ~m)) * Subst(Antiderivative((~x) ^ ~m * (~a * (~x) ^ ~j + ~b * (~x) ^ ~n) ^ ~p, ~x), ~x, ~v) <-- FreeQ([~a, ~b, ~j, ~m, ~n, ~p], ~x) && LinearPairQ(~u, ~v, ~x)
end

max_stars_repo_path: src/rules/1 Algebraic functions/1.1 Binomial products/1.1.4 Improper/1.1.4.2 (c x)^m (a x^j+b x^n)^p.jl


text:
module Keras

include("read.jl")
include("graph.jl")
include("functional.jl")
include("support_functions.jl")

end


text:
using Compat

function gdbm_open(name::String, flag::String="r")
  handle = ccall((:gdbm_open, "libgdbm"), Ptr{Cvoid}, (Cstring, Int32, Int32, Int32, Ptr{Cvoid}), name, 0, OPENFLAGS[flag], 420, C_NULL)
  handle == C_NULL && error("File could not be opened!")
  handle
end

function gdbm_close(handle::Ptr{Cvoid})
  ccall((:gdbm_close, "libgdbm"), Cvoid, (Ptr{Cvoid},), handle)
end

function gdbm_store(handle::Ptr{Cvoid}, key::Datum, value::Datum, flag::String="r")
  ret = ccall((:gdbm_store, "libgdbm"), Int32, (Ptr{Cvoid}, Datum, Datum, Int32), handle, key, value, STOREFLAGS[flag])
  ret == -1 && error("Database is not writable or key/value is not a valid string.")
  ret == 1 && error("Key is already in database.")
  nothing
end

function gdbm_fetch(handle::Ptr{Cvoid}, key::Datum)
  datum = ccall((:gdbm_fetch, "libgdbm"), Datum, (Ptr{Cvoid}, Datum), handle, key)
  datum.dptr == C_NULL && throw(KeyError(key))
  datum
end

function gdbm_exists(handle::Ptr{Cvoid}, key::Datum)
  ret = ccall((:gdbm_exists, "libgdbm"), Int32, (Ptr{Cvoid}, Datum), handle, key)
  ret == 0 && return false
  true
end

function gdbm_count(handle::Ptr{Cvoid})
  count = Ref(UInt(0))
  ret = ccall((:gdbm_count, "libgdbm"), Int32, (Ptr{Cvoid}, Ref{UInt}), handle, count)
  ret == -1 && error("Error reading database.")
  Int(count[])
end

function gdbm_delete(handle::Ptr{Cvoid}, key::Datum)
  ret = ccall((:gdbm_delete, "libgdbm"), Int32, (Ptr{Cvoid}, Datum), handle, key)
  ret ≠ 0 && error("Database is not writable or key not found.")
  nothing
end

function gdbm_firstkey(handle::Ptr{Cvoid})
  ccall((:gdbm_firstkey, "libgdbm"), Datum, (Ptr{Cvoid}, ), handle)
end

function gdbm_nextkey(handle::Ptr{Cvoid}, prev::Datum)
  ccall((:gdbm_nextkey, "libgdbm"), Datum, (Ptr{Cvoid}, Datum), handle, prev)
end

function libc_free(ptr::Ptr{UInt8})
  ccall((:free, "libc"), Cvoid, (Ptr{UInt8}, ), ptr)
end


text:
################################################################################
# Douglas-Rachford splitting iterator
#
# See:
# [1] Eckstein, Bertsekas "On the Douglas-Rachford Splitting Method and the Proximal Point Algorithm for Maximal Monotone Operators*", Mathematical Programming, vol. 55, no. 1, pp. 293-318 (1989).
#

struct DRSIterator{I <: Integer, R <: Real, T <: AbstractArray{R}} <: ProximalAlgorithm{I,T}
    x::T
    f
    g
    gamma::R
    maxit::I
    tol::R
    verbose::I
    verbose_freq::I
    y::T
    r::T
    z::T
    FPR_x::T
end

################################################################################
# Constructor(s)

function DRSIterator(x0::AbstractArray{R}; f=Zero(), g=Zero(), gamma::R=1.0, maxit::I=10000, tol::R=1e-4, verbose=1, verbose_freq = 100) where {I, R}
    y = copy(x0)
    r = copy(x0)
    z = copy(x0)
    FPR_x = copy(x0)
    FPR_x .= Inf
    return DRSIterator{I, R, typeof(x0)}(x0, f, g, gamma, maxit, tol, verbose, verbose_freq, y, r, z, FPR_x)
end

################################################################################
# Utility methods

maxit(sol::DRSIterator) = sol.maxit

converged(sol::DRSIterator, it) = it > 0 && maximum(abs,sol.FPR_x)/sol.gamma <= sol.tol

verbose(sol::DRSIterator) = sol.verbose > 0
verbose(sol::DRSIterator, it) = sol.verbose > 0 && (sol.verbose == 2 ? true : (it == 1 || it%sol.verbose_freq == 0))

function display(sol::DRSIterator)
    @printf("%6s | %10s | %10s |\n ", "it", "gamma", "fpr")
    @printf("------|------------|------------|\n")
end

function display(sol::DRSIterator, it)
    @printf("%6d | %7.4e | %7.4e |\n", it, sol.gamma, maximum(abs,sol.FPR_x)/sol.gamma)
end

function Base.show(io::IO, sol::DRSIterator)
    println(io, "Douglas-Rachford Splitting" )
    println(io, "fpr        : $(maximum(abs,sol.FPR_x))")
    print(  io, "gamma      : $(sol.gamma)")
end

################################################################################
# Initialization

function initialize!(sol::DRSIterator)
    return sol.z
end

################################################################################
# Iteration

function iterate!(sol::DRSIterator{I, T}, it::I) where {I, T}
    prox!(sol.y, sol.f, sol.x, sol.gamma)
    sol.r .= 2 .*sol.y .- sol.x
    prox!(sol.z, sol.g, sol.r, sol.gamma)
    sol.FPR_x .= sol.y .- sol.z
    sol.x .-= sol.FPR_x
    return sol.z
end

################################################################################
# Solver interface

function DRS(x0; kwargs...)
    sol = DRSIterator(x0; kwargs...)
    it, point = run!(sol)
    return (it, point, sol)
end


text:
using Test
#include("../algoritmos-busqueda/generar-vecino.jl")
#include("../algoritmos-busqueda/busqueda-local.jl")
include("../../learner/BL-euclidean-1-NN.jl")

@testset "Buscar vecino cambiando todo " begin
    w = [0.5, 0.5, 0.5]
    m = 0
    s = 0.3
    w_n = GenNeighbourhood(w, m, s)
    for i in 1:3
        @test size(w_n) == size(w)
        @test !(false in map( x-> 0<= x <= 1, w_n))
        w = w_n
        GenNeighbourhood(w, m, s)
    end
end
@testset "Buscar vecino un índice concreto" begin
    for i in 1:5
        w = [0.5, 0.5, 0.5, 0.5, 0.5]
        m = 0
        s = 0.3
        w_n = GenNeighbourhoodIesimo(w, m, s, i)
        @test size(w_n) == size(w)
        @test !(false in map( x-> 0<= x <= 1, w_n))
        @test w[1:i-1] == w_n[1:i-1]
        @test w[i+1:5] == w_n[i+1:5]
    end
end

@testset "Primero mejor" begin
    F(x)=sum(x) # Función de evaluación simple

    for tam in 12:15
        w = PrimeroMejor(10, 2, tam, 0.3, F)
        @test length(w) == tam
        @test !(false in map( x-> 0<= x <= 1, w))
    end
end

@testset "Busqueda local learner" begin
    # El último parámetro es "aleatorio"
    data = [-1 0 1; 1 1 1; -1 -3 1; 1 -2 2; 9 -1 2; 1 2 2; 10 20 1]
    class = ['-', '+',  '-',   '-',  '-', '+' , '+' ]
    clasificador, funcion_evaluacion , w = BL_LearnerOneNN(data, class)

end


text:
n_scen = 10_000       # number of Monte Carlo scenarios
costs_fixed = 20      # fixed costs
s = 0.03              # risk free rate of interest
α = 0.99              # confidence level of economic capital
seed = 2              # random seed for repeatable calculations

insurance_input =
  DataFrame(
  id =          [:fire, :liab, :theft],
  name =        ["Fire", "Liability", "Theft"],
  premium =     [600.,  300.,  100.],
  loss_ratio =  [0.75,  0.75,  0.75],
  cost_ratio =  [0.05,  0.05,  0.05],
  var_coeff =   [0.50,  0.60,  0.70],
  re_ceded =    [0.25,  0.20,  0.20],
  re_costs =    [-0.06, -0.06, -0.06])
## add a column with an explicit counter: the lines of business
insurance_input[:ctr] = collect(1:nrow(insurance_input))

invest_input =
  DataFrame(
  id =    [:invest],
  name = ["Investment"],
  init = [1400.],
  cost_ratio = [0.005],
  mean = [0.05],
  sd =[0.02])
## add a column with an explicit counter: investments
invest_input[:ctr] = [nrow(insurance_input)+1]

tau_kendall = Real[1.0  0.3  0.2  0.0;
                   0.3  1.0  0.6  0.0;
                   0.2  0.6  1.0  0.0;
                   0.0  0.0  0.0  1.0]

hurdle = 0.1          # hurdle rate for EVA calculation
n_cloud = 100         # number of points in cloud diagrams


text:
#=
Some of the below test functions have been modified from their original form to to allow for
tunable input/output sizes, or to test certain programmatic behaviors. Thus, one should not
expect these functions to be "correct" for their original purpose.
=#

########################
# f(x::Number)::Number #
########################

num2num_1(x) = sin(x)^2 / cos(x)^2
num2num_2(x) = 2*x + sqrt(x*x*x)
num2num_3(x) = 10.31^(x + x) - x
num2num_4(x) = x
num2num_5(x) = 1
num2num_6(x) = 1. / (1. + exp(-x))

#######################
# f(x::Number)::Array #
#######################

function num2arr_1(x)
    return reshape([num2num_1(x),
                    num2num_2(x),
                    num2num_3(x),
                    num2num_1(x) - num2num_2(x),
                    num2num_2(x),
                    num2num_3(x),
                    num2num_2(x),
                    num2num_3(x)], 2, 2, 2)
end

########################
# f(x::Vector)::Number #
########################

vec2num_1(x) = (exp(x[1]) + log(x[3]) * x[4]) / x[5]
vec2num_2(x) = x[1]*x[2] + sin(x[1])
vec2num_3(x) = vecnorm(x' .* x)
vec2num_4(x) = ((sum(x) + prod(x)); 1)
vec2num_5(x) = sum((-x).^3)

function rosenbrock_1(x)
    a = one(eltype(x))
    b = 100 * a
    result = zero(eltype(x))
    for i in 1:length(x)-1
        result += (a - x[i])^2 + b*(x[i+1] - x[i]^2)^2
    end
    return result
end

function rosenbrock_2(x)
    a = x[1]
    b = 100 * a
    v = map((i, j) -> (a - j)^2 + b*(i - j^2)^2, x[2:end], x[1:end-1])
    return sum(v)
end

rosenbrock_3(x) = sum(map((i, j) -> (1 - j)^2 + 100*(i - j^2)^2, x[2:end], x[1:end-1]))

function rosenbrock_4(x)
    t1 = (1 + x[1:end-1]).^2
    t2 = x[2:end] + (x[1:end-1]).^2
    return sum(t1 + 100 * (@compat(abs.(t2)).^2))
end

function ackley(x)
    a, b, c = 20.0, -0.2, 2.0*π
    len_recip = inv(length(x))
    sum_sqrs = zero(eltype(x))
    sum_cos = sum_sqrs
    for i in x
        sum_cos += cos(c*i)
        sum_sqrs += i^2
    end
    return (-a * exp(b * sqrt(len_recip*sum_sqrs)) -
            exp(len_recip*sum_cos) + a + e)
end

self_weighted_logit(x) = inv(1.0 + exp(-dot(x, x)))

########################
# f(x::Matrix)::Number #
########################

mat2num_1(x) = det(first(x) * inv(x * x) + x)

function mat2num_2(x)
    a = reshape(x, length(x), 1)
    b = reshape(copy(x), 1, length(x))
    return trace(@compat(log.((1 .+ (a * b)) .+ a .- b)))
end

function mat2num_3(x)
    k = length(x)
    N = isqrt(k)
    A = reshape(x, N, N)
    return sum(map(n -> sqrt(abs(n) + n^2) * 0.5, A))
end

mat2num_4(x) = mean(sum(@compat(sin.(x)) * x, 2))

softmax(x) = sum(@compat(exp.(x)) ./ sum(@compat(exp.(x)), 2))

###########################################
# f(::Matrix, ::Matrix, ::Matrix)::Number #
###########################################

relu(x) = @compat(log.(1.0 .+ @compat(exp.(x))))
sigmoid(n) = 1. / (1. + @compat(exp.(-n)))
neural_step(x1, w1, w2) = sigmoid(dot(w2[1:size(w1, 2)], relu(w1 * x1[1:size(w1, 2)])))

################################
# f!(y::Array, x::Array)::Void #
################################

# Credit for `chebyquad!`, `brown_almost_linear!`, and `trigonometric!` goes to
# Kristoffer Carlsson (@KristofferC).

function chebyquad!(y, x)
    tk = 1/length(x)
    for j = 1:length(x)
        temp1 = 1.0
        temp2 = 2x[j]-1
        temp = 2temp2
        for i = 1:length(y)
            y[i] += temp2
            ti = temp*temp2 - temp1
            temp1 = temp2
            temp2 = ti
        end
    end
    iev = -1.0
    for k = 1:length(y)
        y[k] *= tk
        if iev > 0
            y[k] += 1/(k^2-1)
        end
        iev = -iev
    end
    return nothing
end

function brown_almost_linear!(y, x)
    c = sum(x) - (length(x) + 1)
    for i = 1:(length(x)-1), j = 1:(length(y)-1)
        y[j] += x[i] + c
    end
    y[length(y)] = prod(x) - 1
    return nothing
end

function trigonometric!(y, x)
    for i in x
        for j in eachindex(y)
            y[j] = cos(i)
        end
    end
    c = sum(y)
    n = length(x)
    for i in x
        for j in eachindex(y)
            y[j] = sin(i) * y[j] + n - c
        end
    end
    return nothing
end

function mutation_test_1!(y, x)
    y[1] = x[1]
    y[1] = y[1] * x[2]
    y[2] = y[2] * x[3]
    y[3] = sum(y)
    return nothing
end

function mutation_test_2!(y, x)
    y[1] *= x[1]
    y[2] *= x[1]
    y[1] *= x[2]
    y[2] *= x[2]
    return nothing
end

######################
# f(x::Array)::Array #
######################

chebyquad(x) = (y = zeros(x); chebyquad!(y, x); return y)

brown_almost_linear(x) = (y = zeros(x); brown_almost_linear!(y, x); return y)

trigonometric(x) = (y = ones(x); trigonometric!(y, x); return y)

mutation_test_1(x) = (y = zeros(x); mutation_test_1!(y, x); return y)

mutation_test_2(x) = (y = ones(x); mutation_test_2!(y, x); return y)

arr2arr_1(x) = (sum(x .* x); zeros(x))

arr2arr_2(x) = x[1, :] .+ x[1, :] .+ first(x)


text:
module MaxHelpingHandCustomSandwichLavaSettingsTrigger

using ..Ahorn, Maple

@mapdef Trigger "MaxHelpingHand/CustomSandwichLavaSettingsTrigger" CustomSandwichLavaSettingsTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight,
    onlyOnce::Bool=false, direction::String="CoreModeBased", speed::Number=20.0)

const directions = String["AlwaysUp", "AlwaysDown", "CoreModeBased"]

const placements = Ahorn.PlacementDict(
    "Custom Sandwich Lava Settings (max480's Helping Hand)" => Ahorn.EntityPlacement(
        CustomSandwichLavaSettingsTrigger,
        "rectangle",
    ),
)

Ahorn.editingOptions(entity::CustomSandwichLavaSettingsTrigger) = Dict{String, Any}(
    "direction" => directions
)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand