text:
"""
    emit!
Emit Julia expression.
"""
function emit! end

emit!(dag::ExprDAG, node::ExprNode, options::Dict; args...) = dag

############################### Function ###############################

function _get_func_name(cursor, options)
    library_name = options["library_name"]
    haskey(options, "library_names") || return library_name
    file_name = get_filename(cursor) |> normpath
    ks = filter(x->endswith(file_name, Regex(x)), collect(keys(options["library_names"])))
    if !isempty(ks)
        library_name = options["library_names"][ks[1]] # FIXME: same file name in different folders?
    end
    return library_name
end

"Get argument names and types. Return (names, types)."
function _get_func_arg(cursor, options, dag)
    # argument names
    conflict_syms = get(options, "function_argument_conflict_symbols", [])
    func_args = get_function_args(cursor)
    # handle unnamed args
    arg_names = Vector{Symbol}(undef, length(func_args))
    for (i, arg) in enumerate(func_args)
        ns = Symbol(name(arg))
        safe_name = make_name_safe(ns)
        safe_name = isempty(safe_name) ? "arg$i" : safe_name
        # handle name collisions
        if haskey(dag.tags, ns) || haskey(dag.ids, ns) || haskey(dag.ids_extra, ns)
            safe_name *= "_"
        elseif safe_name ∈ conflict_syms
            safe_name = "_" * safe_name
        end
        arg_names[i] = Symbol(safe_name)
    end

    # argument types
    arg_types = [getArgType(getCursorType(cursor), i - 1) for i in 1:length(func_args)]
    args = Union{Expr,Symbol}[translate(tojulia(arg), options) for arg in arg_types]
    for (i, arg) in enumerate(args)
        # array function arguments should decay to pointers
        # e.g. double f[3] => Ptr{Cdouble} instead of NTuple{3, Cdouble}
        if Meta.isexpr(arg, :curly) && first(arg.args) == :NTuple
            args[i] = Expr(:curly, :Ptr, last(arg.args))
        end

        # do the same test for the underlying type of a typedef
        c = getTypeDeclaration(arg_types[i])
        if c isa CLTypedefDecl
            ty = getCanonicalType(getTypedefDeclUnderlyingType(c))
            jl = translate(tojulia(ty), options)
            if Meta.isexpr(jl, :curly) && first(jl.args) == :NTuple
                args[i] = Expr(:curly, :Ptr, last(jl.args))
            end
        end
    end
    return arg_names, args
end

function _get_func_return_type(cursor, options)
    translate(tojulia(getCursorResultType(cursor)), options)
end

function emit!(dag::ExprDAG, node::ExprNode{FunctionProto}, options::Dict; args...)
    cursor = node.cursor

    library_name = _get_func_name(cursor, options)
    library_expr = Meta.parse(library_name)

    func_name = Symbol(spelling(cursor))

    arg_names, args = _get_func_arg(cursor, options, dag)
    ret_type = _get_func_return_type(cursor, options)

    is_strict_typed = get(options, "is_function_strictly_typed", false)
    signature = if is_strict_typed
        efunsig(func_name, arg_names, args)
    else
        Expr(:call, func_name, arg_names...)
    end

    use_ccall_macro = get(options, "use_ccall_macro", false)
    if use_ccall_macro
        name_type_pairs = [Expr(:(::), n, t) for (n, t) in zip(arg_names, args)]
        body = Expr(:macrocall, Symbol("@ccall"), nothing,
                    Expr(:(::), Expr(:call,
                                        :($library_expr.$func_name),
                                        name_type_pairs...,
                                    ),
                         ret_type
                        )
                    )
    else
        body = :(ccall(
            ($(QuoteNode(func_name)), $library_expr),
            $ret_type,
            $(Expr(:tuple, args...)),
            $(arg_names...),
        ))
    end

    push!(node.exprs, Expr(:function, signature, Expr(:block, body)))

    return dag
end

# TODO: clean up this
function is_ptr_type_expr(@nospecialize t)
    (t === :Cstring || t === :Cwstring) && return true
    isa(t, Expr) || return false
    t = t::Expr
    return t.head === :curly && t.args[1] === :Ptr
end

function efunsig(name::Symbol, args::Vector{Symbol}, types)
    x = [is_ptr_type_expr(t) ? a : Expr(:(::), a, t) for (a, t) in zip(args, types)]
    return Expr(:call, name, x...)
end

function emit!(dag::ExprDAG, node::ExprNode{FunctionNoProto}, options::Dict; args...)
    use_ccall_macro = get(options, "use_ccall_macro", false)
    cursor = node.cursor
    @warn "No Prototype for $cursor - assuming no arguments"
    @assert getNumArguments(cursor) == 0

    library_name = _get_func_name(cursor, options)
    library_expr = Meta.parse(library_name)

    func_name = Symbol(spelling(cursor))
    ret_type = translate(tojulia(getCursorResultType(cursor)), options)

    signature = Expr(:call, func_name)
    if use_ccall_macro
        call = :($library_expr.$func_name()::$ret_type)
        body = Expr(:macrocall, Symbol("@ccall"), nothing, call)
    else
        body = :(ccall(($(QuoteNode(func_name)), $library_expr), $ret_type, $(Expr(:tuple))))
    end

    push!(node.exprs, Expr(:function, signature, Expr(:block, body)))

    return dag
end

function emit!(dag::ExprDAG, node::ExprNode{FunctionVariadic}, options::Dict; args...)
    # @ccall is needed to support variadic argument
    use_ccall_macro = get(options, "use_ccall_macro", true)
    wrap_variadic_function = get(options, "wrap_variadic_function", false)
    if use_ccall_macro && wrap_variadic_function
        cursor = node.cursor
        is_strict_typed = get(options, "is_function_strictly_typed", false)
        arg_names, args = _get_func_arg(cursor, options, dag)
        ret_type = _get_func_return_type(cursor, options)
        library_name = _get_func_name(cursor, options)
        library_expr = Meta.parse(library_name)
        func_name = Symbol(spelling(cursor))

        signature = is_strict_typed ? efunsig(func_name, arg_names, args) : Expr(:call, func_name, arg_names...)
        push!(signature.args, :(va_list...))

        fixed_args = map(arg_names, args) do name, type
            :($name::$type)
        end
        va_list_expr = Expr(:$, :(to_c_type_pairs(va_list)...))
        ccall_body = :($library_expr.$func_name($(fixed_args...); $va_list_expr)::$ret_type)
        ccall_expr = Expr(:macrocall, Symbol("@ccall"), nothing, ccall_body)
        body = quote
            $(Meta.quot(ccall_expr))
        end

        generated = Expr(:function, signature, body)
        ex = Expr(:macrocall, Symbol("@generated"), nothing, generated)
        rm_line_num_node!(ex)
        push!(node.exprs, ex)
    end
    return dag
end

############################### Typedef ###############################

function emit!(dag::ExprDAG, node::ExprNode{TypedefElaborated}, options::Dict; args...)
    if haskey(dag.tags, node.id) || haskey(dag.ids_extra, node.id)
        # generate nothing for duplicated typedef nodes
        # pass
    else
        ty = getTypedefDeclUnderlyingType(node.cursor)
        typedefee = translate(tojulia(ty), options)
        typedef_sym = make_symbol_safe(node.id)
        push!(node.exprs, :(const $typedef_sym = $typedefee))
    end
    return dag
end

function replace_pointee!(expr, s::Symbol)
    if Meta.isexpr(expr, :curly)
        if expr.args[2] isa Expr
            replace_pointee!(expr.args[2], s)
        elseif expr.args[2] isa Symbol
            expr.args[2] = s
        end
    end
    return nothing
end

function emit!(dag::ExprDAG, node::ExprNode{TypedefMutualRef}, options::Dict; args...)
    if haskey(dag.tags, node.id) || haskey(dag.ids_extra, node.id)
        # generate nothing for duplicated typedef nodes
        # pass
    else
        ty = getTypedefDeclUnderlyingType(node.cursor)
        jlty = tojulia(ty)
        leaf_ty = get_jl_leaf_type(jlty)

        # first, we generate a fake struct
        real_sym = leaf_ty.sym
        fake_sym = Symbol("__JL_$real_sym")
        push!(node.exprs, Expr(:struct, true, fake_sym, Expr(:block)))
        # then, generate overloading methods to emulates the behavior
        # `Base.unsafe_load`
        signature = Expr(:call, :(Base.unsafe_load), :(x::Ptr{$fake_sym}))
        body = Expr(:block, :(unsafe_load(Ptr{$real_sym}(x))))
        push!(node.exprs, Expr(:function, signature, body))
        # `Base.getproperty`
        signature = Expr(:call, :(Base.getproperty), :(x::Ptr{$fake_sym}), :(f::Symbol))
        body = Expr(:block, :(getproperty(Ptr{$real_sym}(x), f)))
        push!(node.exprs, Expr(:function, signature, body))
        # `Base.setproperty!`
        signature = Expr(:call, :(Base.setproperty!), :(x::Ptr{$fake_sym}), :(f::Symbol), :v)
        body = Expr(:block, :(setproperty!(Ptr{$real_sym}(x), f, v)))
        push!(node.exprs, Expr(:function, signature, body))

        # generate typedef
        typedefee = translate(jlty, options)
        replace_pointee!(typedefee, fake_sym)
        typedef_sym = make_symbol_safe(node.id)
        push!(node.exprs, :(const $typedef_sym = $typedefee))
    end
    return dag
end

function emit!(dag::ExprDAG, node::ExprNode{TypedefFunction}, options::Dict; args...)
    ty = getTypedefDeclUnderlyingType(node.cursor)
    typedefee = translate(tojulia(ty), options)
    typedef_sym = make_symbol_safe(node.id)
    push!(node.exprs, :(const $typedef_sym = $typedefee))
    return dag
end

function emit!(dag::ExprDAG, node::ExprNode{TypedefDefault}, options::Dict; args...)
    ty = getTypedefDeclUnderlyingType(node.cursor)
    typedefee = translate(tojulia(ty), options)
    typedef_sym = make_symbol_safe(node.id)
    push!(node.exprs, :(const $typedef_sym = $typedefee))
    return dag
end

function emit!(dag::ExprDAG, node::ExprNode{TypedefToAnonymous}, options::Dict; args...)
    ty = getTypedefDeclUnderlyingType(node.cursor)
    typedefee = translate(tojulia(ty), options)
    typedef_sym = make_symbol_safe(node.id)
    push!(node.exprs, :(const $typedef_sym = $typedefee))
    return dag
end

############################### Struct ###############################

function _emit_getproperty_ptr!(body, root_cursor, cursor, options)
    field_cursors = fields(getCursorType(cursor))
    field_cursors = isempty(field_cursors) ? children(cursor) : field_cursors
    for field_cursor in field_cursors
        n = name(field_cursor)
        if isempty(n)
            _emit_getproperty_ptr!(body, root_cursor, field_cursor, options)
            continue
        end
        fsym = make_symbol_safe(n)
        fty = getCursorType(field_cursor)
        ty = translate(tojulia(fty), options)
        offset = getOffsetOf(getCursorType(root_cursor), n)
        if isBitField(field_cursor)
            w = getFieldDeclBitWidth(field_cursor)
            @assert w <= 32 # Bit fields should not be larger than int(32 bits)
            d, r = divrem(offset, 32)
            ex = :(f === $(QuoteNode(fsym)) && return (Ptr{$ty}(x + $(4d)), $r, $w))
        else
            d = offset ÷ 8
            ex = :(f === $(QuoteNode(fsym)) && return Ptr{$ty}(x + $d))
        end
        push!(body.args, ex)
    end
end

# Base.getproperty(x::Ptr, f::Symbol) -> Ptr
function emit_getproperty_ptr!(dag, node, options)
    sym = make_symbol_safe(node.id)
    signature = Expr(:call, :(Base.getproperty), :(x::Ptr{$sym}), :(f::Symbol))
    body = Expr(:block)
    _emit_getproperty_ptr!(body, node.cursor, node.cursor, options)
    push!(body.args, :(return getfield(x, f)))
    getproperty_expr = Expr(:function, signature, body)
    push!(node.exprs, getproperty_expr)
    return dag
end

function rm_line_num_node!(ex::Expr)
    filter!(ex.args) do arg
        arg isa Expr && rm_line_num_node!(arg)
        !isa(arg, LineNumberNode)
    end
    return ex
end

function emit_getproperty!(dag, node, options)
    sym = make_symbol_safe(node.id)

    ref_expr = :(r = Ref{$sym}(x))
    conv_expr = :(ptr = Base.unsafe_convert(Ptr{$sym}, r))
    fptr_expr = :(fptr = getproperty(ptr, f))

    load_expr = :(GC.@preserve r unsafe_load(fptr))
    load_expr.args[2] = nothing

    load_base_expr = :(GC.@preserve r unsafe_load(baseptr32))
    load_base_expr.args[2] = nothing
    load_next_expr = :(GC.@preserve r unsafe_load(baseptr32 + 4))
    load_next_expr.args[2] = nothing

    if is_bitfield_type(node.type)
        ex = quote
            if fptr isa Ptr
                return $load_expr
            else
                baseptr, offset, width = fptr
                ty = eltype(baseptr)
                baseptr32 = convert(Ptr{UInt32}, baseptr)
                u64 = $load_base_expr
                if offset + width > 32
                    u64 |= ($load_next_expr) << 32
                end
                u64 = (u64 >> offset) & ((1 << width) - 1)
                return u64 % ty
            end
        end
    else
        ex = load_expr
    end

    rm_line_num_node!(ex)

    signature = Expr(:call, :(Base.getproperty), :(x::$sym), :(f::Symbol))
    body = Expr(:block, ref_expr, conv_expr, fptr_expr, ex)
    getproperty_expr = Expr(:function, signature, body)

    push!(node.exprs, getproperty_expr)

    return dag
end

function emit_setproperty!(dag, node, options)
    sym = make_symbol_safe(node.id)
    signature = Expr(:call, :(Base.setproperty!), :(x::Ptr{$sym}), :(f::Symbol), :v)
    store_expr = :(unsafe_store!(getproperty(x, f), v))
    if is_bitfield_type(node.type)
        body = quote
            fptr = getproperty(x, f)
            if fptr isa Ptr
                $store_expr
            else
                baseptr, offset, width = fptr
                baseptr32 = convert(Ptr{UInt32}, baseptr)
                u64 = unsafe_load(baseptr32)
                straddle = offset + width > 32
                if straddle
                    u64 |= unsafe_load(baseptr32 + 4) << 32
                end
                mask = ((1 << width) - 1)
                u64 &= ~(mask << offset)
                u64 |= (unsigned(v) & mask) << offset
                unsafe_store!(baseptr32, u64 & typemax(UInt32))
                if straddle
                    unsafe_store!(baseptr32 + 4, u64 >> 32)
                end
            end
        end
        rm_line_num_node!(body)
    else
        body = Expr(:block, store_expr)
    end
    setproperty_expr = Expr(:function, signature, body)
    push!(node.exprs, setproperty_expr)
    return dag
end

function get_names_types(root_cursor, cursor, options)
    field_cursors = fields(getCursorType(cursor))
    field_cursors = isempty(field_cursors) ? children(cursor) : field_cursors
    tys = []
    fsyms = []
    for field_cursor in field_cursors
        n = name(field_cursor)
        if isempty(n)
            _emit_getproperty_ptr!(root_cursor, field_cursor, options)
            continue
        end
        fsym = make_symbol_safe(n)
        fty = getCursorType(field_cursor)
        ty = translate(tojulia(fty), options)
        push!(tys, ty)
        push!(fsyms, fsym)
    end
    fsyms, tys
end

function emit_constructor!(dag, node::ExprNode{<:AbstractUnionNodeType}, options)
    sym = make_symbol_safe(node.id)
    fsyms, tys = get_names_types(node.cursor, node.cursor, options)
    union_sym = Symbol(:__U_, sym)
    push!(node.exprs, :(const $union_sym = Union{$(tys...)}))
    body = Expr(:block,
        :(ref = Ref{$sym}()),
        :(ptr = Base.unsafe_convert(Ptr{$sym}, ref)),
    )

    if get(options, "union_single_constructor", false)
        branch_args = map(zip(fsyms, tys)) do (fsym, ty)
            cond = :(val isa $ty)
            assign = :(ptr.$fsym = val)
            (cond, assign)
        end

        first_args, rest = Iterators.peel(branch_args)
        ex = Expr(:if, first_args...)
        push!(body.args, ex)
        foreach(rest) do (cond, assign)
            _ex = Expr(:elseif, cond, assign)
            push!(ex.args, _ex)
            ex = _ex
        end

        push!(body.args, :(ref[]))
        push!(node.exprs, Expr(:function, Expr(:call, sym, :(val::$union_sym)), body))
    else
        foreach(zip(fsyms, tys)) do (fsym, ty)
            _body = copy(body)
            append!(_body.args,
                [
                    :(ptr.$fsym = $fsym),
                    :(ref[]),
                ]
            )
            func = Expr(:function, :($sym($fsym::$ty)), _body)
            push!(node.exprs, func)
        end
    end
end

function emit_constructor!(dag, node::ExprNode{<:StructLayout}, options)
    sym = make_symbol_safe(node.id)
    fsyms, tys = get_names_types(node.cursor, node.cursor, options)
    body = quote
        ref = Ref{$sym}()
        ptr = Base.unsafe_convert(Ptr{$sym}, ref)
        $((:(ptr.$fsym = $fsym) for fsym in fsyms)...)
        ref[]
    end

    rm_line_num_node!(body)

    func = Expr(:function, Expr(:call, sym, (:($fsym::$ty) for (fsym, ty) in zip(fsyms, tys))...), body)
    push!(node.exprs, func)
end

function emit!(dag::ExprDAG, node::ExprNode{<:AbstractStructNodeType}, options::Dict; args...)
    struct_sym = make_symbol_safe(node.id)
    block = Expr(:block)
    expr = Expr(:struct, false, struct_sym, block)
    field_cursors = fields(getCursorType(node.cursor))
    field_cursors = isempty(field_cursors) ? children(node.cursor) : field_cursors
    for field_cursor in field_cursors
        field_sym = make_symbol_safe(name(field_cursor))
        field_ty = getCursorType(field_cursor)
        push!(block.args, Expr(:(::), field_sym, translate(tojulia(field_ty), options)))
    end
    push!(node.exprs, expr)

    if startswith(string(node.id), "##Ctag") || startswith(string(node.id), "__JL_Ctag")
        emit_getproperty_ptr!(dag, node, options)
        emit_getproperty!(dag, node, options)
        emit_setproperty!(dag, node, options)
    end

    if haskey(options, "field_access_method_list")
        if string(node.id) in options["field_access_method_list"]
            emit_getproperty_ptr!(dag, node, options)
            # emit_getproperty!(dag, node, options)
            emit_setproperty!(dag, node, options)
        end
    end

    return dag
end

function emit!(dag::ExprDAG, node::ExprNode{StructMutualRef}, options::Dict; args...)
    node_idx = args[:idx]
    struct_sym = make_symbol_safe(node.id)
    block = Expr(:block)
    expr = Expr(:struct, false, struct_sym, block)
    mutual_ref_field_cursors = CLCursor[]
    field_cursors = fields(getCursorType(node.cursor))
    field_cursors = isempty(field_cursors) ? children(node.cursor) : field_cursors
    for field_cursor in field_cursors
        field_sym = make_symbol_safe(name(field_cursor))
        field_ty = getCursorType(field_cursor)
        jlty = tojulia(field_ty)
        leaf_ty = get_jl_leaf_type(jlty)
        translated = translate(jlty, options)

        if jlty != leaf_ty && !is_jl_basic(leaf_ty)
            # this assumes tag-types and identifiers that have the same name are the same
            # thing, which is validated in the audit pass.
            field_idx = get(dag.tags, leaf_ty.sym, typemax(Int))
            if field_idx == typemax(Int)
                field_idx = get(dag.ids, leaf_ty.sym, typemax(Int))
            end

            # if `leaf_ty.sym` can not be found in `tags` and `ids` then it's in `ids_extra`
            field_idx == typemax(Int) && @assert haskey(dag.ids_extra, leaf_ty.sym)

            if node_idx < field_idx
                # this assumes that circular references were removed at pointers
                @assert is_jl_pointer(jlty)

                # also emit the original expressions, so we can add corresponding comments
                # in the pretty-print pass
                comment = Expr(:(::), field_sym, deepcopy(translated))
                replace_pointee!(translated, :Cvoid)
                push!(mutual_ref_field_cursors, field_cursor)
                # Avoid pushing two expressions for one field
                push!(block.args, Expr(:block, comment, :($field_sym::$translated)))
                continue
            end
        end

        push!(block.args, Expr(:(::), field_sym, translated))
    end

    push!(node.exprs, expr)

    # make corrections by overloading `Base.getproperty` for those `Ptr{Cvoid}` fields
    if !isempty(mutual_ref_field_cursors)
        getter = Expr(:call, :(Base.getproperty), :(x::$struct_sym), :(f::Symbol))
        body = Expr(:block)
        for mrfield_cursor in mutual_ref_field_cursors
            n = name(mrfield_cursor)
            @assert !isempty(n)
            fsym = make_symbol_safe(n)
            fty = getCursorType(mrfield_cursor)
            ty = translate(tojulia(fty), options)
            ex = :(f === $(QuoteNode(fsym)) && return $ty(getfield(x, f)))
            push!(body.args, ex)
        end
        push!(body.args, :(return getfield(x, f)))
        getproperty_expr = Expr(:function, getter, body)
        push!(node.exprs, getproperty_expr)
    end

    if haskey(options, "field_access_method_list")
        if string(node.id) in options["field_access_method_list"]
            emit_getproperty_ptr!(dag, node, options)
            # emit_getproperty!(dag, node, options)
            emit_setproperty!(dag, node, options)
        end
    end

    return dag
end

# codegen for memory-pool-like structures
function emit!(dag::ExprDAG, node::ExprNode{<:RecordLayouts}, options::Dict; args...)
    sym = make_symbol_safe(node.id)
    n = getSizeOf(getCursorType(node.cursor))
    expr = Expr(:struct, false, sym, Expr(:block, :(data::NTuple{$n,UInt8})))
    push!(node.exprs, expr)

    # emit Base.getproperty(x::Ptr, f::Symbol) etc.
    emit_getproperty_ptr!(dag, node, options)
    emit_getproperty!(dag, node, options)
    emit_setproperty!(dag, node, options)

    opt = get(options, "add_record_constructors", [])
    if (opt isa Bool && opt) || sym in opt
        emit_constructor!(dag, node, options)
    end

    return dag
end

############################### Enum ###############################
const ENUM_SYMBOL2TYPE = Dict(
    :Cchar => Cchar,
    :Cuchar => Cuchar,
    :Cshort => Cshort,
    :Cushort => Cushort,
    :Clong => Clong,
    :Culong => Culong,
    :Cwchar_t => Cwchar_t,
    :Clonglong => Clonglong,
    :Culonglong => Culonglong,
    :Cint => Cint,
    :Cuint => Cuint,
    :Cintmax_t => Cintmax_t,
    :Cuintmax_t => Cuintmax_t,
    :Cptrdiff_t => Cptrdiff_t,
    :Cssize_t => Cssize_t,
    :UInt64 => UInt64,
    :UInt32 => UInt32,
    :UInt16 => UInt16,
    :UInt8 => UInt8,
    :Int64 => Int64,
    :Int32 => Int32,
    :Int16 => Int16,
    :Int8 => Int8,
)

function emit!(dag::ExprDAG, node::ExprNode{<:AbstractEnumNodeType}, options::Dict; args...)
    cursor = node.cursor
    ty = translate(tojulia(getEnumDeclIntegerType(cursor)), options)
    # there is no need to treat __attribute__ specially as we directly query the integer
    # type of the enum from Clang
    enum_type = ENUM_SYMBOL2TYPE[ty]
    @assert Core.sizeof(enum_type) == getSizeOf(getCursorType(cursor))

    enum_sym = make_symbol_safe(node.id)
    expr = Expr(:macrocall, Symbol("@cenum"), nothing, Expr(:(::), enum_sym, enum_type))
    push!(node.exprs, expr)

    for item_cursor in children(cursor)
        item_cursor isa CLEnumConstantDecl || continue
        item_name = spelling(item_cursor)
        @assert !isempty(item_name)
        item_sym = make_symbol_safe(item_name)
        val = value(item_cursor)
        push!(node.exprs, :($item_sym = $val))
    end

    return dag
end

################## ForwardDecl OpaqueDecl DefaultDecl DuplicatedTags ##################

# skip non-opaque forward decls
emit!(dag::ExprDAG, node::ExprNode{<:ForwardDecls}, options::Dict; args...) = dag

function emit!(dag::ExprDAG, node::ExprNode{<:OpaqueTags}, options::Dict; args...)
    struct_sym = make_symbol_safe(node.id)
    if get(options, "opaque_as_mutable_struct", true)
        push!(node.exprs, Expr(:struct, true, struct_sym, Expr(:block)))
    else
        push!(node.exprs, :(const $struct_sym = Cvoid))
    end
    return dag
end

function emit!(dag::ExprDAG, node::ExprNode{<:UnknownDefaults}, options::Dict; args...)
    @error "missing implementation for $(node), please file an issue to Clang.jl."
    return dag
end

# skip duplicated nodes
emit!(dag::ExprDAG, node::ExprNode{<:DuplicatedTags}, options::Dict; args...) = dag


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

JLLWrappers.@generate_main_file_header("libwatcher")
JLLWrappers.@generate_main_file("libwatcher", UUID("f2e3f4ed-6d30-53f6-bb33-6a6ab454fa9c"))
end  # module libwatcher_jll

ext: jl


text:
using Pkg
path = @__DIR__
cd(path)
Pkg.activate(path)

using Mitosis
using MitosisStochasticDiffEq
using LinearAlgebra, Statistics, Random, StatsBase
using DelimitedFiles
using StaticArrays
using NewickTree
using DifferentialEquations
using DiffEqNoiseProcess

const d = 2
const 𝕏 = SVector{d,Float64}
const MS = Mitosis
const MSDE = MitosisStochasticDiffEq

# using MeasureTheory
# import MeasureTheory.logdensity
MS.dim(p::WGaussian{(:μ,:Σ,:c)}) = length(p.μ)
MS._logdet(p::WGaussian{(:μ,:Σ,:c)}) = MS._logdet(p.Σ, MS.dim(p))
MS.whiten(p::WGaussian{(:μ,:Σ,:c)}, x) = MS.lchol(p.Σ)\(x - p.μ)
MS.sqmahal(p::WGaussian, x) = MS.norm_sqr(MS.whiten(p, x))
MS.logdensity(p::WGaussian{(:μ,:Σ,:c)}, x) = p.c - (MS.sqmahal(p,x) + MS._logdet(p) + MS.dim(p)*log(2pi))/2

include("tree.jl")
include("sdetree.jl")

Random.seed!(10)

## Read tree
tree = [Tree(S50), Tree(S20coal), Tree(S)][1]
#print_tree(tree.newick)

## define model (two traits, as in presentation)
bθ = @SVector [.5, 1.2]
σ0 = @SVector [log(0.2), log(1.05)]
θ0 = (bθ, σ0)

#const M = [-1.0 1.0; 1.0 -1.0]
const M = SMatrix{d,d}(-1.0, 1.0, 1.0, -1.0)
f(u,θ,t) = 𝕏(tanh.(Diagonal(θ[1]) * M * u) )  # f(u,θ,t) = Diagonal(θ[1]) * M * u
g(u,θ,t) = Diagonal(exp.(θ[2]))

## forward sample on the tree to simulate the observations at the tips
dt0 = 0.001
u0 = zero(𝕏)  # value at root node
Xd, segs = forwardsample(tree, u0, θ0, dt0, f, g, EM(false))

## Define auxiliary process
B(θ) = Diagonal(θ[1]) * M
##Σ(θ) = Diagonal(exp.(θ[2]))
#σ̃(θ) = SMatrix{d,d}(exp(θ[2][1]), 0.0, 0.0, exp(θ[2][2])) # or?  g(0,θ,0) * g(0, θ, 0)'
σ̃(θ) = g(0,θ,0)

function mcmc2(tree, Xd, f, g, θinit, prior;
                 ρ=0.99,
                 iters=5000,
                 dt=0.01,
                 σprop=0.05,
                 precisionatleaves=10e5,
                 apply_time_change=true,
                 𝒫=(:μ,:Σ,:c),  # 𝒫=(:F,:\Gamma,:c)
                 recomputeguidingterm=true,
                 alg=Tsit5(),
                 SDEalg=EM(false)
                 )
    #σ0 = θinit[2]

    if 𝒫==(:μ,:Σ,:c)
        #Q = [i in tree.lids ? WGaussian{(:μ,:Σ,:c)}(Vector(Xd[i]), precisionatleaves*Matrix(I(d)), 0.0) : missing for i in tree.ids]
        # in case of static arrays
        leavescov = inv(precisionatleaves) * SA_F64[1 0; 0 1]
        Q = [i in tree.lids ? WGaussian{(:μ,:Σ,:c)}(Xd[i], leavescov, 0.0) : missing for i in tree.ids]
    elseif 𝒫==(:F,:Γ,:c)
        leavesprecision = precisionatleaves * SA_F64[1 0; 0 1]
        Q = [i in tree.lids ? WGaussian{(:F,:Γ,:c)}(leavesprecision*Xd[i], leavesprecision, 0.0) : missing for i in tree.ids]
    else
        @error "𝒫 not defined"
    end
    θ = θinit
    θlin = (B(θinit), zeros(d), σ̃(θinit))

    Q, messages = bwfiltertree!(Q, tree, θlin, dt; apply_time_change=apply_time_change, alg=alg)
    Qᵒ, messagesᵒ = Q, messages

    guidedsegs = Vector{Any}(undef, tree.n) # save all guided segments
    X = zeros(𝕏, tree.n)  # values at nodes
    for id in tree.lids
        X[id] = Xd[id]
    end

    Z = [myinnov(messages[i].ts, 𝕏) for i ∈ 2:tree.n]

    X, guidedsegs, ll, 𝐋 = fwguidtree!(X, guidedsegs, Q, messages, tree, f, g, θ, Z, SDEalg)

    Xᵒ, guidedsegsᵒ, Qᵒ = deepcopy(X), deepcopy(guidedsegs), deepcopy(Q)
    θs = [θinit]
    accepted = 0

    for iter in 1:iters
        #θᵒ = (θ[1] + σprop * randn(size(θ[1])...), σ0)
        θᵒ = (θ[1] + σprop * randn(size(θ[1])...), θ[2] + σprop * randn(size(θ[2])...))
        Zᵒ =  [pcn_innov(Z[i], ρ, 𝕏) for i ∈ eachindex(Z)]
        θlinᵒ = (0.9B(θᵒ), zeros(d), σ̃(θᵒ))
        if recomputeguidingterm==true
            Qᵒ, messagesᵒ = bwfiltertree!(Qᵒ, tree, θlinᵒ, dt, apply_time_change=apply_time_change,  alg=alg)
        else
            Qᵒ = Q
        end
        Xᵒ, guidedsegsᵒ, llᵒ, 𝐋ᵒ = fwguidtree!(Xᵒ, guidedsegsᵒ, Qᵒ, messagesᵒ, tree, f, g, θᵒ, Zᵒ, SDEalg)
        Δ = 𝐋ᵒ - 𝐋 + logdensity(prior[1], θᵒ[1]) - logdensity(prior[1], θ[1]) + logdensity(prior[2], θᵒ[2]) - logdensity(prior[2], θ[2])

        if mod(iter, 10) == 0
             println(iter,"   ", round.(θ[1];digits=2),"   ", round.(θᵒ[1];digits=2), "       ",
                    round(𝐋;digits=3), "   ", round(𝐋ᵒ;digits=3),"   ", round(Δ; digits=3))
        end

        if log(rand()) < Δ
            𝐋 = 𝐋ᵒ
            ll .= llᵒ
            θ = θᵒ
            θlin = θlinᵒ
            guidedsegs, guidedsegsᵒ = guidedsegsᵒ, guidedsegs # don't care about this
            X, Xᵒ = Xᵒ, X # don't care about this
            Z, Zᵒ = Zᵒ, Z    #  Z .= Zᵒ
            Q, Qᵒ = Qᵒ, Q
            accepted += 1
        end
        push!(θs, deepcopy(θ))
    end

    θs, guidedsegs, accepted/iters, (Xᵒ, guidedsegsᵒ, Qᵒ, messagesᵒ, tree)
end

prior = (MS.Gaussian{(:F,:Γ)}(zeros(2), Matrix(0.01*I(2))) , MS.Gaussian{(:F,:Γ)}(-ones(2), Matrix(0.01*I(2))))
#θinit = (SVector(3.2, 1.0), σ0)
θinit = (SVector(3.2, 1.0), SVector(.1, .1))

iters = 50_000
@time θs, guidedsegs, frac_accepted, forwardguiding_input = mcmc2(tree, Xd, f, g, θinit, prior;
  iters=iters)#, 𝒫=(:F,:Γ,:c))#, dt = dt0)

## summary stats

println("fraction accepted equals: $frac_accepted")
burnin = div(3iters,4):iters
θs1 = getindex.(θs,1)[burnin]
@show mean(getindex.(θs1,1))
@show mean(getindex.(θs1,2))
@show θ0[1][1]
@show θ0[1][2]
println("-------")
θs2 = getindex.(θs,2)[burnin]
@show mean(getindex.(θs2,1))
@show mean(getindex.(θs2,2))
@show θ0[2][1]
@show θ0[2][2]

## plotting
PLOT = true
if PLOT==true
    include("plottingtree.jl")
end

## partially conjugate steps
STOP = true
if STOP==false
    ## regression of drift parameters θ[1]
    function driftparamstree(prior, segs, tree::Tree, f, g, θ, paramjac, messages)
        G = deepcopy(prior) # posterior that should be returned
        for i in eachindex(tree.T)
            i == 1 && continue  # skip root-node  (has no parent)
            ts = messages[i].ts
            κ = MSDE.SDEKernel(f, g, ts, θ)
            R = MSDE.Regression(κ,θ=θ[1],paramjac=paramjac)
            G = MSDE.conjugate(R, map(x->x[1:end-1], segs[i].u), G, ts)  # map because last element is the loglikelihood (this is is we provide segs as guidedsegs)
        end
        G
    end

    # Gaussian prior
    prior = MS.Gaussian{(:F,:Γ)}(zeros(2), Matrix(0.001*I(2)))
    # paramjac for non AD version
    function paramjac(u,p,t)
      Diagonal(M*u)
    end
    G = driftparamstree(prior, guidedsegs, tree, f, g, θ, paramjac, messages)
    p̂ = mean(G)
    se = sqrt.(diag(cov(G)))
    display(map((p̂, se, p) -> "$(round(p̂, digits=3)) ± $(round(se, digits=3)) (true: $p)", p̂, se, θ0[1]))
end

# function mcmc(tree, Xd, f, g, θlin, θinit, prior; ρ=0.99, iters=5000, dt=0.01, σprop=0.05, precisionatleaves=10e-6)
#     σ0 = θinit[2]
#     Q = [i in tree.lids ? WGaussian{(:μ,:Σ,:c)}(Vector(Xd[i]), precisionatleaves*Matrix(I(d)), 0.0) : missing for i in tree.ids]
#     Q, messages = bwfiltertree!(Q, tree, θlin, dt)
#
#     guidedsegs = Vector{Any}(undef, tree.n) # save all guided segments
#     X = zeros(𝕏, tree.n)  # values at nodes
#     for id in tree.lids
#         X[id] = Xd[id]
#     end
#
#     Z = [innov(messages[i].ts) for i ∈ 2:tree.n]
#     θ = θinit
#     X, guidedsegs, ll, 𝐋 = fwguidtree!(X, guidedsegs, Q, messages, tree, f, g, θ, Z)
#
#     Xᵒ, guidedsegsᵒ = deepcopy(X), deepcopy(guidedsegs)
#     θs = [θinit]
#     accepted = 0
#
#     for iter in 1:iters
#         θᵒ = (θ[1] + σprop * randn(size(θ[1])...), σ0)
#         Zᵒ = [pcn_innov(Z[i], ρ) for i ∈ eachindex(Z)]
#         Xᵒ, guidedsegsᵒ, llᵒ, 𝐋ᵒ = fwguidtree!(Xᵒ, guidedsegsᵒ, Q, messages, tree, f, g, θᵒ, Zᵒ)
#         Δ = 𝐋ᵒ - 𝐋 + logdensity(prior, θᵒ[1]) - logdensity(prior, θ[1])
#
#         if mod(iter, 100) == 0
#              println(iter,"   ", round.(θ[1];digits=2),"   ", round.(θᵒ[1];digits=2), "       ",
#                     round(𝐋;digits=3), "   ", round(𝐋ᵒ;digits=3),"   ", round(Δ; digits=3))
#         end
#
#         if log(rand()) < Δ
#             𝐋 = 𝐋ᵒ
#             ll .= llᵒ
#             θ = θᵒ
#             guidedsegs, guidedsegsᵒ = guidedsegsᵒ, guidedsegs # don't care about this
#             X, Xᵒ = Xᵒ, X # don't care about this
#             Z, Zᵒ = Zᵒ, Z    #  Z .= Zᵒ
#             accepted += 1
#         end
#         push!(θs, deepcopy(θ))
#     end
#
#     θs, guidedsegs, accepted/iters
# end


text:
# ---
# title: 453. Minimum Moves to Equal Array Elements
# id: problem453
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Math
# link: <https://leetcode.com/problems/minimum-moves-to-equal-array-elements/description/>
# hidden: true
# ---
#
# Given a **non-empty** integer array of size _n_ , find the minimum number of
# moves required to make all array elements equal, where a move is incrementing
# _n_ \- 1 elements by 1.
#
# **Example:**
#
#
#
#     Input:
#     [1,2,3]
#
#     Output:
#     3
#
#     Explanation:
#     Only three moves are needed (remember each move increments two elements):
#
#     [1,2,3]  =>  [2,3,3]  =>  [3,4,3]  =>  [4,4,4]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using GalerkinSparseGrids
using Test
using HCubature
using StaticArrays
using ODE
using SparseArrays

#--------------------------------------
# Testing PDE solvers
#--------------------------------------

# The wave equation, testing conservation of energy

# In the position basis:

import GalerkinSparseGrids.energy_func_1D
import GalerkinSparseGrids.energy_func

@testset "solvers.jl" begin
    @info "Testing wave equation solver 1-D position DG basis... "

    global k = 4
    global level = 4
    global f0 = x->sin(2*pi*x[1])
    global v0 = x->2*pi*cos(2*pi*x[1])
    pos_soln = wave_evolve_1D(k, level, f0, v0, 0, 1; basis="pos", order="45") #problem is here
    energy_soln = energy_func_1D(k, level, pos_soln; basis="pos")
    @testset "Wave equation solver 1D position DG basis" begin
        for i in 1:length(energy_soln[2])
            @test sqrt(energy_soln[2][i]) ≈ 2*pi atol=1.0e-7
        end
    end

    # In the hierarchical basis:
    @info "Testing wave equation solver 1-D hierarchical DG basis... "
    @testset "Wave equation solver 1D hier DG basis" begin
        hier_soln  = wave_evolve_1D(k, level, f0, v0, 0, 1; basis="hier", order="45")
        energy_soln = energy_func_1D(k, level, hier_soln; basis="hier")
        for i in 1:length(energy_soln[2])
            @test sqrt(energy_soln[2][i]) ≈ 2*pi atol=1.0e-7
        end

        hier_soln  = wave_evolve(1, k, level, f0, v0, 0, 1; order="78")
        energy_soln = energy_func(1, k, level, hier_soln)
        for i in 1:length(energy_soln[2])
            @test sqrt(energy_soln[2][i]) ≈ 2*pi atol=1.0e-7
        end
    end

    # 2D Sparse case:

    @info "Testing wave equation solver on 2-D sparse DG basis... "

    @testset "Wave equation solver 2D sparse DG basis" begin

        D = 2
        k_used = 3
        n_used = 5

        f0 = x->sin(2*pi*x[1])*sin(2*pi*x[2])
        v0 = x->0
        sparse_soln = wave_evolve(D, k_used, n_used, f0, v0, 0, 1; order="78", scheme="sparse");
        senergy = energy_func(D, k_used, n_used, sparse_soln)
        @test senergy[2][1]-senergy[2][end] > 0
        @test abs(senergy[2][1]-senergy[2][end]) < 1.0e-8
        for i in 1:length(senergy[1])
            @test sqrt(senergy[2][i]) ≈ sqrt(2)*pi atol=1.0e-4
        end

        sparse_soln = wave_evolve(D, k_used, n_used, f0, v0, 0,1; order="45", scheme="sparse");
        senergy = energy_func(D, k_used, n_used, sparse_soln)
        @test senergy[2][1] > senergy[2][end]
        @test senergy[2][1] ≈ senergy[2][end] atol=1.0e-8
        for i in 1:length(senergy[1])
            @test sqrt(senergy[2][i]) ≈ sqrt(2)*pi atol=1.0e-4
        end
    end

    # @info "Testing traveling wave example in 2-D... "
    # @testset "Traveling wave example 2D" begin
    #     m = [1, 2]
    #     truesoln = x -> cos(2*pi*(dot(m,x) - sqrt(dot(m,m))*0.54))
    #
    #     D = 2
    #     k_used = 3
    #     n_used = 6
    #
    #     soln = traveling_wave_solver(k_used, n_used, m, 0, 0.54; order="45")
    #     dict = V2D(D, k_used, n_used, soln[2][end]; scheme="sparse")
    #     @test mcerr(x->reconstruct_DG(dict, [x...]), truesoln, D) < 0.05
    #
    #     soln = traveling_wave_solver(k_used, n_used, m, 0, 0.54; order="78")
    #     dict = V2D(D, k_used, n_used, soln[2][end]; scheme="sparse")
    #     @test mcerr(x->reconstruct_DG(dict, [x...]), truesoln, D) < 0.05
    # end
end


text:
# an example on mixed integer programming
#
#   maximize x + 2 y + 5 z
#
#   s.t.  x + y + z <= 10
#         x + 2 y + z <= 15
#
#         x is continuous: 0 <= x <= 5
#         y is integer: 0 <= y <= 10
#         z is binary
#

using CPLEX
using JuMP

m = Model(solver=CplexSolver())
@defVar(m, 5 >= x >= 0)
@defVar(m, 10 >= y >= 0, Int)
@defVar(m, z, Bin)
@setObjective(m, Max, x + 2*y + 5*z)
@addConstraint(m, x + y + z <= 10)
@addConstraint(m, x + 2*y + z <= 15)
status = solve(m)
println("walltime = $(MathProgBase.getsolvetime(getInternalModel(m)))")


text:
"""
   Plan

Stores core elements of the state of the current districting plan: the assignment
of each node, the population of each district, and the set of cut edges (the edges
of the plan's graph that define the district boundaries). Intended to be updated
whenever a step of a chain is accepted.
"""
mutable struct Plan
    n_districts::Int
    assignment::Array{Int}
    district_populations::Array{Int}
    district_nodes::Array{BitSet}
    district_edges::Array{BitSet}  # per-district edge lists
    district_adj::Array{Int, 2}    # cut edge count (seam length) between districts
    cut_edges::BitSet
    n_cut_edges::Int
end

function Plan(graph::IndexedGraph, assignment::Array{Int})::Plan
    n_districts = maximum(assignment) - minimum(assignment) + 1
    district_populations = zeros(Int, n_districts)
    for (index, node_pop) in enumerate(graph.population)
        district_populations[assignment[index]] += node_pop
    end
    cut_edges = BitSet(Int[])
    district_adj = zeros(Int,  n_districts, n_districts)
    district_edges = [BitSet() for _ in 1:n_districts]
    district_nodes = [BitSet() for _ in 1:n_districts]
    for index in 1:graph.n_edges
        left_assignment = assignment[graph.edges[1, index]]
        right_assignment = assignment[graph.edges[2, index]]
        if left_assignment != right_assignment
            push!(cut_edges, index)
            district_adj[left_assignment, right_assignment] += 1
            district_adj[right_assignment, left_assignment] += 1
        end
        push!(district_edges[left_assignment], index)
        push!(district_edges[right_assignment], index)
        push!(district_nodes[left_assignment], graph.edges[1, index])
        push!(district_nodes[right_assignment], graph.edges[2, index])
    end
    return Plan(n_districts, assignment, district_populations,
                district_nodes, district_edges, district_adj,
                cut_edges, length(cut_edges))
end

"""
    Plan(graph, assignment_col)

Create a new plan using data stored in `graph`. The assignment of each node
is loaded from the column `assignment_col` in the metadata of `graph`.
"""
function Plan(graph::IndexedGraph, assignment_col::AbstractString)::Plan
    assignment = zeros(Int, graph.n_nodes)
    for index in 1:graph.n_nodes
        assignment[index] = graph.attributes[index][assignment_col]
    end
    return Plan(graph, assignment)
end


text:
function makeearlieststartday(N, K, C, S)
    ret = zeros(Int, K)
    d = 1
    k = 1
    while k <= K
        if S[d] == 'o'
            ret[k] = d
            d += C + 1
            k += 1
        elseif S[d] == 'x'
            d += 1
        end
    end
    ret
end

function makelateststartday(N, K, C, S)
    ret = zeros(Int, K)
    d = N
    k = K
    while k >= 1
        if S[d] == 'o'
            ret[k] = d
            d -= C + 1
            k -= 1
        else S[d] == 'x'
            d -= 1
        end
    end
    ret
end

function solve()
    N, K, C = [parse(Int, x) for x in split(readline())]
    S = readline()

    earlieststartday = makeearlieststartday(N, K, C, S)
    lateststartday = makelateststartday(N, K, C, S)
    for (e, l) in zip(earlieststartday, lateststartday)
        if e == l
            println(e)
        end
    end
end

solve()

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
function InnerProduct(in1,in2)

	return convert(Float32,real(sum(conj(in1[:]).*in2[:])))

end

function InnerProduct(in1::String,in2::String;ntrace=10000)
	# Inner Product of two vectors with same length

	#	ip = 0.0
	#	itrace_in = 1
	#	nx = GetNumTraces(in1)
	#	while itrace_in <= nx
	#		d1,h1,e = SeisRead(in1,group="some",itrace=itrace_in,ntrace=ntrace)
	#		d2,h2,e = SeisRead(in2,group="some",itrace=itrace_in,ntrace=ntrace)
	#		ip += InnerProduct(d1,d2)
	#		num_traces_in = size(d1,2)
	#		itrace_in += num_traces_in
	#	end

	#	return ip

	d1,h1,ext = SeisRead(in1)
	d2,h2,ext = SeisRead(in2)

	return sum(d1[:].*d2[:])

end

function InnerProduct(in1::Array{String,1},in2::Array{String,1};ntrace=10000)

		ip = 0.0
		for j = 1 : length(in1)
			ip += InnerProduct(in1[j],in2[j],ntrace=ntrace)
		end

		return ip

end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 2/20/2015

# Multiple aircraft encounter model based on generating encounters pairwise using CorrAEM
# DBN samples are generated at runtime at each step.

module PairwiseCorrAEMDBNImpl

export
    AddObserver,

    getInitialState,
    initialize,
    update,
    get,

    PairwiseCorrAEMDBN

using AbstractEncounterDBNImpl
using AbstractEncounterDBNInterfaces
using CommonInterfaces
using Util
using Encounter
using RLESUtils, Observers
using Base.Test

import CommonInterfaces.addObserver
import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractEncounterDBNInterfaces.get
import AbstractEncounterDBNInterfaces.getInitialState

import CorrAEMImpl
AEM_CorrAEMParameters = CorrAEMImpl.CorrAEMParameters
AEM_getInitialState = CorrAEMImpl.getInitialState
AEM_generateEncounter = CorrAEMImpl.generateEncounter

include(Pkg.dir("SISLES/src/Encounter/CorrAEMImpl/corr_aem_sample.jl"))

type PairwiseCorrAEMDBN <: AbstractEncounterDBN

  number_of_aircraft::Int64

  parameter_file::ASCIIString
  parameters::AEM_CorrAEMParameters

  dirichlet_initial
  dirichlet_transition

  encounter_seed::UInt64

  aem::CorrAEM
  ac1_initial_dist::Union{Vector{Union{Function, Void}}, Void}
  ac1_dynamic_vars::Union{Array{Float64,2},Void}

  t::Float64

  #initial state
  initial_states::Vector{CorrAEMInitialState}

  #initial command
  initial_commands_d::Vector{Vector{Int64}} #discrete current state
  initial_commands::Vector{Vector{Float64}} #continuous of current dynamic variables

  #current command
  commands_d::Vector{Vector{Int64}} #discrete current state
  commands::Vector{Vector{Float64}} #continuous of current dynamic variables

  #pre-allocated output to avoid repeating reallocations
  output_commands::Vector{CorrAEMCommand}

  logProb::Float64 #log probability of output

  function PairwiseCorrAEMDBN(number_of_aircraft::Int, parameter_file::AbstractString,
                              encounter_seed::UInt64;
                              initial_file::AbstractString="",
                              transition_file::AbstractString="")

    dbn = new()

    @test number_of_aircraft >= 2 #multi-aircraft sim
    dbn.number_of_aircraft     = number_of_aircraft

    dbn.parameter_file = parameter_file
    dbn.encounter_seed = encounter_seed

    const default_number_of_initial_samples = 1
    const default_number_of_transition_samples = 50

    #Based on CorrAEM (call it repeatedly for pairwise encounters)
    dbn.aem = CorrAEM(parameter_file, initial_file, default_number_of_initial_samples,
                      transition_file, default_number_of_transition_samples,
                      b_read_from_file=false, b_write_to_file=false)

    dbn.parameters = dbn.aem.parameters
    dbn.dirichlet_initial = bn_dirichlet_prior(dbn.parameters.N_initial)
    dbn.dirichlet_transition = bn_dirichlet_prior(dbn.parameters.N_transition)

    dbn.ac1_initial_dist = nothing
    dbn.ac1_dynamic_vars = nothing

    dbn.t = 0.0

    srand(encounter_seed) #There's a rand inside generateEncounter
    generate_encounter(dbn) #sets initial_states, initial_commands_d, initial_commands

    dbn.output_commands = CorrAEMCommand[ CorrAEMCommand(0.,0.,0.,0.) for i = 1:number_of_aircraft ]
    dbn.logProb = 0.0

    return dbn
  end
end

const AC1_MASK = [true,true,false,false,true,false,true,false,true,false,true,false,true,false,false,false]
const SAVE_INDICES_AC1 = [1, 10, 12, 14]

const map_G2L = Dict(2 => 1, 9 => 2, 11 => 3, 13 => 4, 17 => 5, 19 => 6) #global to local
const map_L2G = Dict(1 => 2, 2 => 9, 3 => 11, 4 => 13, 5=> 17, 6 => 19) #local to global
const map_var2ind_L = Dict(:L => 1, :v_d => 2, :h_d0 => 3, :psi_d0 => 4, :h_d1 => 5, :psi_d1 => 6) #variable names to local
const map_ind2var_L = Dict(1 => :L, 2 => :v_d, 3 => :h_d0, 4 => :psi_d0, 5 => :h_d1, 6 => :psi_d1) #local to variable names
const temporal_map = [3 5; 4 6] #[dynamic_variables0; dynamic_variables1]

function generate_encounter(dbn::PairwiseCorrAEMDBN; sample_number::Int64 = 0)

  @test dbn.number_of_aircraft >= 2

  p = dbn.parameters
  dynamic_variables1 = temporal_map[:,2]

  #initial aircraft states and commands
  dbn.initial_states = Array(CorrAEMInitialState,dbn.number_of_aircraft)
  dbn.initial_commands = Array(Vector{Float64},dbn.number_of_aircraft) #[L,v_d,h_d,psi_d]
  dbn.initial_commands_d = Array(Vector{Int64},dbn.number_of_aircraft) #[L,v_d,h_d,psi_d,[hd_tp1,psid_tp1]]

  #Handle the first 2 aircraft
  AEM_generateEncounter(dbn.aem)

  for i = 1:2 #aircraft_number
    dbn.initial_states[i] = AEM_getInitialState(dbn.aem,i)

    t, v_d, h_d, psi_d = getNextCommand(dbn.aem, i)
    dbn.initial_commands[i] = Float64[dbn.aem.L, v_d, h_d, psi_d]
    initial_commands_d = discretize(dbn.parameters,unconvert_units(dbn.initial_commands[i]))

    dbn.initial_commands_d[i] = [ initial_commands_d, Int64(zeros(dynamic_variables1)) ]
  end

  #Save aircraft 1 initial vars
  dbn.ac1_initial_dist = masked_initial_dist(dbn.aem.initial,AC1_MASK)

  #Save aircraft 1 dynamic vars
  dbn.ac1_dynamic_vars = save_ac1_dyn_vars(squeeze(dbn.aem.states[1,:,:],1))

  #Handle aircraft 3+
  for i = 3:dbn.number_of_aircraft

    #Load aircraft 1 initial vars into initial distribution
    dbn.aem.initial_distributions = dbn.ac1_initial_dist

    states = get_em_sample(dbn.aem)

    #load aircraft 1 dynamic vars into states (replacing this one)
    load_dynamic_vars!(states,dbn.ac1_dynamic_vars)

    sim_and_transform(dbn.aem,states)

    #aircraft 1 is globally common, aircraft 2 is newly generated and stacked on
    dbn.initial_states[i] = AEM_getInitialState(dbn.aem,2)

    t, v_d, h_d, psi_d = getNextCommand(dbn.aem, 2)
    dbn.initial_commands[i] = Float64[dbn.aem.L, v_d, h_d, psi_d]
    initial_commands_d = discretize(dbn.parameters,unconvert_units(dbn.initial_commands[i]))

    dbn.initial_commands_d[i] = [ initial_commands_d, Int64(zeros(dynamic_variables1)) ]
  end
end

#Excerpted from CorrAEMImpl.jl and made into its own function
function get_em_sample(aem::CorrAEM)

  if aem.initial_distributions == nothing
    states = em_sample(aem, aem.number_of_transition_samples)
  else
    states, aem.ISInfo = em_sample(aem, aem.number_of_transition_samples, initial_dist = aem.initial_distributions)
  end

  aem.number_of_encounters_generated += 1

  states = states[:, 2:end]
  convert_units(states) #moved from sim_and_transform

  return states
end

function load_ac1_dyn_vars!(states::Array{Float64,3},dyn_vars)
  for (l,g) in enumerate(SAVE_INDICES_AC1)
    states[1,:, g] = dyn_vars[:, l]
  end
end

#Excerpted from CorrAEMImpl.jl and made into its own function
function sim_and_transform(aem::CorrAEM,states)
  #convert_units(states) #moved to sim_and_transform

  t, A, L, chi, beta_, C1, C2, v1, v2, v1d, v2d, h1d, h2d, psi1d, psi2d, hmd, vmd = states[1, :]

  aem.initial = vec(states[1, 2:end])

  aem.A = A
  aem.L = L
  aem.C[1] = C1
  aem.C[2] = C2
  aem.v_init[1] = v1
  aem.v_init[2] = v2

  aem.geometry_at_TCA = [chi, beta_, hmd, vmd]

  aem.states[1, :, 1] = states[:, 1]  # time
  aem.states[1, :, 2] = states[:, 10] # v1d
  aem.states[1, :, 3] = states[:, 12] # h1d
  aem.states[1, :, 4] = states[:, 14] # psi1d
  aem.state_index[1] = 0

  aem.states[2, :, 1] = states[:, 1]  # time
  aem.states[2, :, 2] = states[:, 11] # v2d
  aem.states[2, :, 3] = states[:, 13] # h2d
  aem.states[2, :, 4] = states[:, 15] # psi2d
  aem.state_index[2] = 0

  aem.dn_state_index[1] = 0
  aem.dn_state_index[2] = 0

  simulate_tracks(aem)

  transform_regarding_TCA(aem, aem.L, aem.geometry_at_TCA)

  for i = 1:aem.number_of_aircraft
    aem.state_index[i] = 0
  end

  #for i = 1:(aem.number_of_transition_samples + 1)
  #    print(reshape(aem.dynamic_states[1, i, :], 6)')
  #end

  #for i = 1:aem.number_of_transition_samples
  #    print(reshape(aem.states[1, i, :], 4)')
  #end
end

function masked_initial_dist(v::Vector{Float64}, mask::Vector{Bool})
  return [m ? ()->(x,0.0) : nothing for (m,x) = zip(mask,v)]
end

addObserver(dbn::PairwiseCorrAEMDBN, f::Function) = add_observer(aem.observer, f)
addObserver(dbn::PairwiseCorrAEMDBN, tag::AbstractString, f::Function) = add_observer(aem.observer, tag, f)

function initialize(dbn::PairwiseCorrAEMDBN)
  #reset to initial state
  dbn.commands_d = deepcopy(dbn.initial_commands_d)
  dbn.commands = deepcopy(dbn.initial_commands)
  dbn.t = 0.0

end

function getInitialState(dbn::PairwiseCorrAEMDBN, index::Int)
  return dbn.initial_states[index]
end

function update(dbn::PairwiseCorrAEMDBN)
  logProb = 0.0 #to accumulate over each aircraft

  for i = 1:dbn.number_of_aircraft
    logProb += step_dbn(dbn.parameters,dbn.dirichlet_transition,
                        dbn.commands_d[i],dbn.commands[i])

    dbn.output_commands[i].t = dbn.t
    dbn.output_commands[i].v_d = dbn.commands[i][map_var2ind_L[:v_d]]
    dbn.output_commands[i].h_d = dbn.commands[i][map_var2ind_L[:h_d0]]
    dbn.output_commands[i].psi_d = dbn.commands[i][map_var2ind_L[:psi_d0]]
  end

  dbn.t += 1.0

  return logProb
end

function step_dbn(p::AEM_CorrAEMParameters, dirichlet_transition,
                         command_d::Vector{Int64}, command::Vector{Float64})
  logProb = 0.0

  dynamic_variables0 = temporal_map[:,1] #[hdot_1(t), psidot_1(t)]
  dynamic_variables1 = temporal_map[:,2] #[hdot_1(t+1), psidot_1(t+1)]

  for (o,i_L) in enumerate(dynamic_variables1)
    i_G = map_L2G[i_L]
    parents_G = p.G_transition[:, i_G]
    if !isempty(find(parents_G))
      parents_L = Int64[map_G2L[iparents] for iparents in find(parents_G)]
      j_G = asub2ind(p.r_transition[parents_G], command_d[parents_L]')
      weights = p.N_transition[i_G][:, j_G] + dirichlet_transition[i_G][:, j_G]
      weights /= sum(weights)
      command_d[i_L] = select_random(weights)
      logProb += log(weights[command_d[i_L]])
      #Resampling and dediscretizing process
      i0_L = dynamic_variables0[o]
      i0_G = map_L2G[i0_L]
      if (command_d[i_L] != command_d[i0_L]) || #compare to state at last time step, #Different bin, do resample
        (command_d[i_L] == command_d[i0_L] && rand() < p.resample_rates[i0_G]) #Same bin but meets resample rate
        command[i0_L] = dediscretize(command_d[i_L],p.boundaries[i0_G],p.zero_bins[i0_G])
        command[i0_L] = convert_units(command[i0_L],map_ind2var_L[i0_L])
      end
      #Else same bin and does not meet rate, just set equal to previous (no update)
    end
  end

  # update x(t) with x(t+1)
  command_d[dynamic_variables0] = command_d[dynamic_variables1]

  #return
  return logProb
end

function get(dbn::PairwiseCorrAEMDBN, aircraft_number::Int)
  return dbn.output_commands[aircraft_number]
end

function val2ind(boundariesi,ri,value)
  if !isempty(boundariesi)
    index = findfirst(x -> (x > value), boundariesi) - 1

    if index == -1
      index = ri
    end
  else
    index = value
  end
  return index
end

function discretize(p::AEM_CorrAEMParameters,v::Vector{Float64})
  return Int64[ val2ind(p.boundaries[map_L2G[i]],
                        p.r_transition[map_L2G[i]],val)
               for (i,val) in enumerate(v) ]
end

function dediscretize(dval::Int64,boundaries::Vector{Float64},zero_bin::Int64)
  val_min = boundaries[dval]
  val_max = boundaries[dval+1]

  return dval == zero_bin ? 0.0 : val_min +  rand() * (val_max - val_min)
end

function asub2ind(siz, x)
# ASUB2IND Linear index from multiple subscripts.
#   Returns a linear index from multiple subscripts assuming a matrix of a
#   specified size.
#
#   NDX = ASUB2IND(SIZ,X) returns the linear index NDX of the element in a
#   matrix of dimension SIZ associated with subscripts specified in X.

    k = [1, cumprod(siz[1:end-1])]
    ndx = k' * (x' - 1) + 1

    return ndx[1]
end

convert_units(v::Vector{Float64}) = Float64[convert_units(v[i],map_ind2var_L[i]) for i=1:endof(v)]
unconvert_units(v::Vector{Float64}) = Float64[unconvert_units(v[i],map_ind2var_L[i]) for i=1:endof(v)]

function convert_units(x::Float64,var::Symbol)
  if var == :v_d0 || var == :v_d1
    return x * 1.68780
  elseif var == :h_d0 || var == :h_d1
    return x / 60
  else
    return x
  end
end

function unconvert_units(x::Float64,var::Symbol)
  if var == :v_d0 || var == :v_d1
    return x / 1.68780
  elseif var == :h_d0 || var == :h_d1
    return x * 60
  else
    return x
  end
end

end #module

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
####
#### Derivative operators
####

@inline ∂xᶜᵃᵃ(i, j, k, grid, u) = δxᶜᵃᵃ(i, j, k, grid, u) / Δx(i, j, k, grid)
@inline ∂xᶠᵃᵃ(i, j, k, grid, c) = δxᶠᵃᵃ(i, j, k, grid, c) / Δx(i, j, k, grid)

@inline ∂yᵃᶜᵃ(i, j, k, grid, v) = δyᵃᶜᵃ(i, j, k, grid, u) / Δy(i, j, k, grid)
@inline ∂yᵃᶠᵃ(i, j, k, grid, c) = δyᵃᶠᵃ(i, j, k, grid, c) / Δy(i, j, k, grid)

@inline ∂zᵃᵃᶜ(i, j, k, grid, w) = δzᵃᵃᶜ(i, j, k, grid, w) / ΔzF(i, j, k, grid)
@inline ∂zᵃᵃᶠ(i, j, k, grid, c) = δzᵃᵃᶠ(i, j, k, grid, c) / ΔzC(i, j, k, grid)

####
#### Derivative operators acting on functions
####

@inline ∂xᶜᵃᵃ(i, j, k, grid, f::F, args...) where F<:Function = (f(i+1, j, k, grid, args...) - f(i,   j, k, grid, args...)) / Δx(i, j, k, grid)
@inline ∂xᶠᵃᵃ(i, j, k, grid, f::F, args...) where F<:Function = (f(i,   j, k, grid, args...) - f(i-1, j, k, grid, args...)) / Δx(i, j, k, grid)

@inline ∂yᵃᶜᵃ(i, j, k, grid, f::F, args...) where F<:Function = (f(i, j+1, k, grid, args...) - f(i, j,   k, grid, args...)) / Δy(i, j, k, grid)
@inline ∂yᵃᶠᵃ(i, j, k, grid, f::F, args...) where F<:Function = (f(i, j,   k, grid, args...) - f(i, j-1, k, grid, args...)) / Δy(i, j, k, grid)

@inline ∂zᵃᵃᶜ(i, j, k, grid, f::F, args...) where F<:Function = (f(i, j, k+1, grid, args...) - f(i, j, k,   grid, args...)) / ΔzF(i, j, k, grid)
@inline ∂zᵃᵃᶠ(i, j, k, grid, f::F, args...) where F<:Function = (f(i, j, k,   grid, args...) - f(i, j, k-1, grid, args...)) / ΔzC(i, j, k, grid)

####
#### Second derivatives
####

@inline ∂²xᶜᵃᵃ(i, j, k, grid, c) = ∂xᶜᵃᵃ(i, j, k, grid, ∂xᶠᵃᵃ, c)
@inline ∂²xᶠᵃᵃ(i, j, k, grid, u) = ∂xᶠᵃᵃ(i, j, k, grid, ∂xᶜᵃᵃ, u)

@inline ∂²yᵃᶜᵃ(i, j, k, grid, c) = ∂yᵃᶜᵃ(i, j, k, grid, ∂yᵃᶠᵃ, c)
@inline ∂²yᵃᶠᵃ(i, j, k, grid, v) = ∂yᵃᶠᵃ(i, j, k, grid, ∂yᵃᶜᵃ, v)

@inline ∂²zᵃᵃᶜ(i, j, k, grid, c) = ∂zᵃᵃᶜ(i, j, k, grid, ∂zᵃᵃᶠ, c)
@inline ∂²zᵃᵃᶠ(i, j, k, grid, w) = ∂zᵃᵃᶠ(i, j, k, grid, ∂zᵃᵃᶜ, w)


text:
module IntroAlgoCh2

function insertionsort!(A::AbstractArray)
    for j=2:length(A)
        key = A[j]

        i = j - 1
        while i > 0 && A[i] > key
            A[i+1] = A[i]
            i = i - 1
        end
        A[i+1] = key
    end
    return A
end

function _merge!(A::Array, p::Int, q::Int, r::Int)
    n1 = q - p + 1
    n2 = r - q

    L = A[p:q]
    R = A[q+1:r]

    i=1
    j=1

    for k=p:r
        if L[i] <= R[j]
            A[k] = L[i]
            i+=1
            if i >  n1
                A[k+1:r] .= R[j:end]
                break
            end
        else
            A[k] = R[j]
            j+=1

            if j > n2
                A[k+1:r] .= L[i:end]
                break
            end
        end
    end
end

function _mergesort!(A::Array, p::Int, r::Int)
    if p<r
        q = (p+r)÷2
        _mergesort!(A,p,q)
        _mergesort!(A,q+1,r)
        _merge!(A,p,q,r)
    end
end

function mergesort!(A::Array)
    _mergesort!(A, 1, length(A))

    return A
end

function mergesort_iterative!(A::Array)
    N = length(A)
    B = similar(A)

    #sorting
    for k=1:2:N-1
        if A[k] <= A[k+1]
            B[k] = A[k]
            B[k+1] = A[k+1]
        else
            B[k] = A[k+1]
            B[k+1] = A[k]
        end
    end
    if isodd(N)
        B[N] = A[N]
    end
    A, B = B, A #A is nu gesorteerd (deels) en B is garbage

    #merging
    mergestep = 2
    while mergestep < N
        for start = 1:mergestep*2:N
            if mergestep >= N - start
                break #er blijft nog maar een partitie over en die is al gesorteerd
            end

            p = start  #eerst index rij1
            q = start + mergestep #eerste index rij2
            r = min(q + mergestep - 1, N) #laatste index rij2

            i = p
            j = q

            for k=p:r
                if A[i] <= A[j]
                    B[k] = A[i]
                    i+=1
                    if i >= q
                        @views B[k+1:r] = A[j:r]
                        break
                    end
                else
                    B[k] = A[j]
                    j+=1

                    if j > r
                        @views B[k+1:r] =  A[i:q-1]
                        break
                    end
                end
            end
            A, B = B, A
        end
        mergestep *= 2
    end

    return A
end

function mergesort_iterative_optimized!(A::Array)
    N = length(A)
    if N <= 64
        return insertionsort!(A)
    end

    B = similar(A)

    #sorting
    sortrange = 1:64:N-64
    for i=sortrange
        insertionsort!(@views A[i:i+64])
    end
    insertionsort!(@views A[maximum(sortrange)+1:N])
    #A is nu gesorteerd (deels) en B is garbage

    #merging
    mergestep = 64
    while mergestep < N
        for start = 1:mergestep*2:N
            if mergestep >= N - start
                break #er blijft nog maar een partitie over en die is al gesorteerd
            end

            p = start  #eerst index rij1
            q = start + mergestep #eerste index rij2
            r = min(q + mergestep - 1, N) #laatste index rij2

            i = p
            j = q

            for k=p:r
                if A[i] <= A[j]
                    B[k] = A[i]
                    i+=1
                    if i >= q
                        @views B[k+1:r] = A[j:r]
                        break
                    end
                else
                    B[k] = A[j]
                    j+=1

                    if j > r
                        @views B[k+1:r] =  A[i:q-1]
                        break
                    end
                end
            end
            A, B = B, A
        end
        mergestep *= 2
    end

    return A
end

function _binarysearch(S_sorted::Vector{Int}, target::Int, low::Int, high::Int)::Union{Nothing, Int}
    if low == high
        return (S_sorted[low] == target ? low : nothing)
    end

    mid = (low + high)÷2

    if S_sorted[mid] == target
        return mid
    elseif S_sorted[mid] > target
        return _binarysearch(S_sorted, target, low, mid)
    else
        return _binarysearch(S_sorted, target, mid, high)
    end
end

function sumfinder(S::Vector{Int}, x::Int)::Union{Nothing, Tuple{Int, Int}}
    mergesort!(S)

    n = length(S)

    for i=1:n-1
        idx = let
            target = x - S[i]
            low = i+1
            high = n

            _binarysearch(S, target, low, high)
        end
        if idx !== nothing
            return (i, idx)
        end
    end

    return nothing

end

export insertionsort!, mergesort!, _binarysearch, sumfinder, mergesort_iterative!, mergesort_iterative_optimized!

end


text:
###
  Name     : c6_13_Pkg_add.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 3/1/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

Pkg.add("AbstractTable")


text:
export Interaction,
    InteractionData,
    ApplicationCommand,
    ApplicationCommandOption,
    ApplicationCommandChoice,
    Component,
    SelectOption,
    ResolvedData

struct ResolvedData
    users::Optional{Dict{Snowflake, User}}
    members::Optional{Dict{Snowflake, Member}}
    roles::Optional{Dict{Snowflake, Role}}
    channels::Optional{Dict{Snowflake, Channel}}
    messages::Optional{Dict{Snowflake, Message}}
end
@boilerplate ResolvedData :constructors :lower :merge

"""
Application Command Choice.
More details [here](https://discord.com/developers/docs/interactions/application-commands#application-command-object-application-command-option-choice-structure).
"""
struct ApplicationCommandChoice
    name::String
    value::Union{String, Number}
end
@boilerplate ApplicationCommandChoice :constructors :docs :lower :merge

# TODO: Custom type gen for `value` field of ApplicationCommandOption.
"""
Application Command Option.
More details [here](https://discord.com/developers/docs/interactions/application-commands#application-command-object-application-command-option-structure).
"""
struct ApplicationCommandOption
    value::Any
    type::Optional{Int}
    name::Optional{String}
    description::Optional{String}
    required::Optional{Bool}
    min_value::Optional{Number}
    max_value::Optional{Number}
    autocomplete::Optional{Bool}
    choices::Optional{Vector{ApplicationCommandChoice}}
    options::Optional{Vector{ApplicationCommandOption}}
    channel_types::Optional{Vector{Int}}
    focused::Optional{Bool}
end
@boilerplate ApplicationCommandOption :constructors :docs :lower :merge

"""
An Application Command.
More details [here](https://discord.com/developers/docs/interactions/application-commands#application-commands).
"""
struct ApplicationCommand <: DiscordObject
    id::OptionalNullable{Snowflake}
    type::Optional{Int}
    application_id::Snowflake
    guild_id::Optional{Snowflake}
    name::String
    description::String
    options::Optional{Vector{ApplicationCommandOption}}
    default_permissions::Optional{Bool}
    version::Optional{Snowflake}
end
@boilerplate ApplicationCommand :constructors :docs :lower :merge :mock

"""
A select option.
More details [here](https://discord.com/developers/docs/interactions/message-components#select-menu-object-select-option-structure).
"""
struct SelectOption <: DiscordObject
    label::String
    value::String
    description::Optional{String}
    emoji::Optional{Emoji}
    default::Optional{Bool}
end
@boilerplate SelectOption :constructors :docs :lower :merge :mock

function SelectOption(label, value, args...)
    if length(args)>0
        r = [:description, :emoji, :default]
        t = Dict([(r[x+1], args[x+1]) for x in 0:length(args)])
        SelectOption(; label=label, value=value, t...)
    else SelectOption(; label=label, value=value) end
end

"""
An interactable component.
More details [here](https://discord.com/developers/docs/interactions/message-components).
"""
struct Component <: DiscordObject
    type::Int
    custom_id::Optional{String}
    value::Optional{String}
    disabled::Optional{Bool}
    style::Optional{Int}
    label::Optional{String}
    emoji::Optional{Emoji}
    url::Optional{String}
    options::Optional{Vector{SelectOption}}
    placeholder::Optional{String}
    min_values::Optional{Int}
    max_values::Optional{Int}
    components::Optional{Vector{Component}}
end
@boilerplate Component :constructors :docs :lower :merge :mock

"""
Data for an interaction.
More details [here](https://discord.com/developers/docs/interactions/receiving-and-responding#interaction-object-interaction-data-structure).
"""
struct InteractionData <: DiscordObject
    id::OptionalNullable{Snowflake}
    name::OptionalNullable{String}
    type::OptionalNullable{Int}
    resolved::Optional{ResolvedData}
    options::Optional{Vector{ApplicationCommandOption}}
    custom_id::OptionalNullable{String}
    component_type::OptionalNullable{Int}
    components::OptionalNullable{Vector{Component}}
    values::Optional{Vector{String}}
    target_id::Optional{Snowflake}
end
@boilerplate InteractionData :constructors :docs :lower :merge

"""
An interaction.
More details [here](https://discord.com/developers/docs/interactions/receiving-and-responding#interaction-object-interaction-structure).
"""
struct Interaction <: DiscordObject
    id::Nullable{Snowflake}
    application_id::Nullable{Snowflake}
    type::Int
    data::OptionalNullable{InteractionData}
    guild_id::Optional{Snowflake}
    channel_id::Optional{Snowflake}
    member::Optional{Member}
    user::Optional{User}
    token::String
    version::Optional{Int}
    message::Optional{Message}
end
@boilerplate Interaction :constructors :docs :lower :merge :mock


text:
module BinaryBuilder
using Libdl, LibGit2, Random, JSON
using BinaryBuilderBase
using ObjectFile
using GitHub
import InteractiveUtils
using Pkg, Base.BinaryPlatforms, Pkg.Artifacts
using ghr_jll

# Re-export useful stuff from Pkg:
export HostPlatform, platform_dlext, valid_dl_path, arch, libc,
       libgfortran_version, libstdcxx_version, cxxstring_abi, detect_libgfortran_version,
       detect_libstdcxx_version, detect_cxxstring_abi, call_abi, wordsize, triplet,
       select_platform, platforms_match, AbstractPlatform, Platform
# Re-export some compat wrappers:
using Pkg.BinaryPlatforms: Linux, MacOS, Windows, FreeBSD
export Linux, MacOS, Windows, FreeBSD
# BinaryBuilderBase/src/Prefix.jl
export Prefix, bindir, libdirs, includedir, logdir, activate, deactivate,
       isinstalled, install, uninstall, list_tarball_files, verify, temp_prefix
# BinaryBuilderBase/src/Rootfs.jl
export supported_platforms, expand_gfortran_versions, expand_cxxstring_abis
# BinaryBuilderBase/src/Platforms.jl
export AnyPlatform
# BinaryBuilderBase/src/Products.jl
export Product, LibraryProduct, FileProduct, ExecutableProduct, FrameworkProduct, satisfied,
       locate, write_deps_file, variable_name
# BinaryBuilderBase/src/Dependency.jl
export Dependency, BuildDependency
# BinaryBuilderBase/src/Sources.jl
export ArchiveSource, FileSource, GitSource, DirectorySource
# Auditor.jl
export audit, collect_files, collapse_symlinks

# Import some useful stuff from BinaryBuilderBase
import BinaryBuilderBase: runshell

include("Auditor.jl")
include("Wizard.jl")

using OutputCollectors, BinaryBuilderBase, .Auditor, .Wizard

include("AutoBuild.jl")
include("Declarative.jl")
include("Logging.jl")

function __init__()
    # If we're running on Azure, enable azure logging:
    if !isempty(get(ENV, "AZP_TOKEN", ""))
        enable_azure_logging()
    end
end

function get_bb_version()
    # Get BinaryBuilder.jl's version and git sha
    version = Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "Project.toml"))["version"]
    try
        # get the gitsha if we can
        repo = LibGit2.GitRepo(dirname(@__DIR__))
        gitsha = string(LibGit2.GitHash(LibGit2.GitCommit(repo, "HEAD")))
        return VersionNumber("$(version)-git-$(gitsha[1:10])")
    catch
        try
            # Settle for the treehash otherwise
            env = Pkg.Types.Context().env
            bb_uuid = Pkg.Types.UUID("12aac903-9f7c-5d81-afc2-d9565ea332ae")
            treehash = bytes2hex(env.manifest[bb_uuid].tree_hash.bytes)
            return VersionNumber("$(version)-tree-$(treehash[1:10])")
        catch
            # Something went so wrong, we can't get any of that.
            return VersionNumber(version)
        end
    end
end

"""
    versioninfo()

Helper function to print out some debugging information
"""
function versioninfo()
    @info("Julia versioninfo(): ")
    InteractiveUtils.versioninfo()

    @info("BinaryBuilder.jl version: $(get_bb_version())")

    @static if Sys.isunix()
        @info("Kernel version: $(readchomp(`uname -r`))")
    end

    # Dump if some important directories are encrypted:
    @static if Sys.islinux()
        print_enc(n, path) = begin
            is_encrypted, mountpoint = is_ecryptfs(path)
            if is_encrypted
                @info("$n is encrypted on mountpoint $mountpoint")
            else
                @info("$n is NOT encrypted on mountpoint $mountpoint")
            end
        end

        print_enc("pkg dir", dirname(@__FILE__))
        print_enc("storage dir", storage_dir())
    end

    # Dump any relevant environment variables:
    @info("Relevant environment variables:")
    env_var_suffixes = [
        "AUTOMATIC_APPLE",
        "USE_SQUASHFS",
        "STORAGE_DIR",
        "RUNNER",
        "ALLOW_ECRYPTFS",
        "USE_CCACHE",
    ]
    for e in env_var_suffixes
        envvar = "BINARYBUILDER_$(e)"
        if haskey(ENV, envvar)
            @info("  $(envvar): \"$(ENV[envvar])\"")
        end
    end

    # Print out the preferred runner stuff here:
    @info("Preferred runner: $(preferred_runner())")

    # Try to run 'echo julia' in Linux x86_64 environment
    @info("Trying to run `echo hello julia` within a Linux x86_64 environment...")

    runner = preferred_runner()(
        pwd();
        cwd="/workspace/",
        platform=Platform("x86_64", "linux"),
        verbose=true
    )
    run_interactive(runner, `/bin/bash -c "echo hello julia"`)

    # If we use ccache, dump the ccache stats
    if BinaryBuilderBase.use_ccache
        @info("ccache stats:")
        runner = preferred_runner()(
            pwd();
            cwd="/workspace/",
            platform=Platform("x86_64", "linux"),
        )
        run_interactive(runner, `/usr/bin/ccache -s`)
    end
    return nothing
end

# Precompilation ahoy!
# include("../deps/SnoopCompile/precompile/precompile_BinaryBuilder.jl")
# _precompile_()

end # module


text:
ArfRef(; prec::Integer = DEFAULT_PRECISION[]) = Arf(; prec = prec)
ArbRef(; prec::Integer = DEFAULT_PRECISION[]) = Arb(; prec = prec)
AcbRef(; prec::Integer = DEFAULT_PRECISION[]) = Acb(; prec = prec)
MagRef() = Mag()

function ArfRef(
    ptr::Ptr{arf_struct},
    parent::Union{arb_struct,ArbRef};
    prec::Integer = DEFAULT_PRECISION[],
)
    ArfRef(ptr, prec, parent)
end
function ArbRef(
    ptr::Ptr{arb_struct},
    parent::Union{acb_struct,AcbRef,arb_vec_struct,arb_mat_struct};
    prec::Integer = DEFAULT_PRECISION[],
)
    ArbRef(ptr, prec, parent)
end
function AcbRef(
    ptr::Ptr{acb_struct},
    parent::Union{Nothing,acb_vec_struct,acb_mat_struct};
    prec::Integer = DEFAULT_PRECISION[],
)
    AcbRef(ptr, prec, parent)
end

Mag(x::MagRef) = set!(Mag(), x)
Arf(x::ArfRef; prec::Integer = precision(x)) = set!(Arf(prec = prec), x)
Arb(x::ArbRef; prec::Integer = precision(x)) = set!(Arb(prec = prec), x)
Acb(x::AcbRef; prec::Integer = precision(x)) = set!(Acb(prec = prec), x)

Base.zero(::Union{Type{MagRef},MagRef}) = zero(Mag)
Base.one(::Union{Type{MagRef},MagRef}) = one(Mag)
for (TRef, T) in [(ArfRef, Arf), (ArbRef, Arb), (AcbRef, Acb)]
    @eval begin
        Base.zero(x::$TRef) = $T(0, prec = precision(x))
        Base.zero(::Type{$TRef}) = zero($T)
        Base.one(x::$TRef) = $T(1, prec = precision(x))
        Base.one(::Type{$TRef}) = one($T)
    end
end

Base.getindex(x::MagRef) = Mag(x)
Base.getindex(x::ArfRef) = Arf(x)
Base.getindex(x::ArbRef) = Arb(x)
Base.getindex(x::AcbRef) = Acb(x)

function midref(x::ArbLike, prec = precision(x))
    mid_ptr = ccall(@libarb(arb_mid_ptr), Ptr{arf_struct}, (Ref{arb_struct},), x)
    ArfRef(mid_ptr, prec, parentstruct(x))
end
function radref(x::ArbLike)
    rad_ptr = ccall(@libarb(arb_rad_ptr), Ptr{mag_struct}, (Ref{arb_struct},), x)
    MagRef(rad_ptr, parentstruct(x))
end

function realref(z::AcbLike; prec = precision(z))
    real_ptr = ccall(@libarb(acb_real_ptr), Ptr{arb_struct}, (Ref{acb_struct},), z)
    ArbRef(real_ptr, prec, parentstruct(z))
end
function imagref(z::AcbLike; prec = precision(z))
    real_ptr = ccall(@libarb(acb_imag_ptr), Ptr{arb_struct}, (Ref{acb_struct},), z)
    ArbRef(real_ptr, prec, parentstruct(z))
end


text:
using PlutoStaticHTMLDemo
using Test

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


text:
module TestVariantCallFormat

using Test
using VariantCallFormat

import BioCore: isfilled
import BioCore.Exceptions.MissingFieldException
import BufferedStreams: BufferedInputStream
import BioCore.Testing: get_bio_fmt_specimens
import YAML

fmtdir = get_bio_fmt_specimens()

include("vcf.jl")
include("bcf.jl")

end


text:
include("../module/CustomStencil.jl")
include("../misc/misc.jl")
include("../misc/cpu_stencils.jl")
## Star Stencil definition with radius = 4
coefs = [1;-0.5;0.25;-0.125;0.05]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 4, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_single = NewStencilInstance(st_def, m_step=false)
st_16bdim =  NewStencilInstance(st_def, m_step=false, bdim=16)

coefs = [1;-0.5;0.25;-0.125;0.05;-0.5;0.25;-0.125;0.05]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 8, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def, m_step=false)
## Input Data size Definition
radius = 4
nx = 9
ny = 9
nz = 9
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

data = CreateData(dx,dy,dz)

t_steps = 4
## GPU Implementation
gpu_out = ApplyStencil(st_inst, data, t_steps)
gpu_out = ApplyStencil(st_single, data, t_steps)
gpu_out = ApplyStencil(st_16bdim, data, t_steps)

CUDA.cuProfilerStart()
gpu_out = ApplyStencil(st_inst, data, t_steps)
gpu_out = ApplyStencil(st_single, data, t_steps)
gpu_out = ApplyStencil(st_16bdim, data, t_steps)

exit()


text:
# ---
# title: 87. Scramble String
# id: problem87
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/scramble-string/description/>
# hidden: true
# ---
#
# We can scramble a string s to get a string t using the following algorithm:
#
#   1. If the length of the string is 1, stop.
#   2. If the length of the string is > 1, do the following:
#     * Split the string into two non-empty substrings at a random index, i.e., if the string is `s`, divide it to `x` and `y` where `s = x + y`.
#     * **Randomly**  decide to swap the two substrings or to keep them in the same order. i.e., after this step, `s` may become `s = x + y` or `s = y + x`.
#     * Apply step 1 recursively on each of the two substrings `x` and `y`.
#
# Given two strings `s1` and `s2` of **the same length** , return `true` if `s2`
# is a scrambled string of `s1`, otherwise, return `false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: s1 = "great", s2 = "rgeat"
#     Output: true
#     Explanation: One possible scenario applied on s1 is:
#     "great" --> "gr/eat" // divide at random index.
#     "gr/eat" --> "gr/eat" // random decision is not to swap the two substrings and keep them in order.
#     "gr/eat" --> "g/r / e/at" // apply the same algorithm recursively on both substrings. divide at ranom index each of them.
#     "g/r / e/at" --> "r/g / e/at" // random decision was to swap the first substring and to keep the second substring in the same order.
#     "r/g / e/at" --> "r/g / e/ a/t" // again apply the algorithm recursively, divide "at" to "a/t".
#     "r/g / e/ a/t" --> "r/g / e/ a/t" // random decision is to keep both substrings in the same order.
#     The algorithm stops now and the result string is "rgeat" which is s2.
#     As there is one possible scenario that led s1 to be scrambled to s2, we return true.
#
#
# **Example 2:**
#
#
#
#     Input: s1 = "abcde", s2 = "caebd"
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: s1 = "a", s2 = "a"
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `s1.length == s2.length`
#   * `1 <= s1.length <= 30`
#   * `s1` and `s2` consist of lower-case English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module TestCompressed
using Accessors: @set
using Test
using PhaseSpaceIO
const PHSP = PhaseSpaceIO
using PhaseSpaceIO.Testing
using PhaseSpaceIO: Format, FormatEGS, FormatIAEA

@testset "==" begin
    for p1 in [arbitrary(IAEAParticle{rand(0:4), rand(0:4)}) for _ in 1:10]
        path = IAEAPath(tempname())
        phsp_write(path, [p1])
        p2 = PhspVector(path)[1]

        @test p2 == p1
        @test p1 == p2
        @test @set(p1.z += 0.1) != p2
        @test p2 != @set p1.E += 0.1
        @test p2 != @set p1.x += 0.1
        @test p2 != @set p1.y += 0.1
        @test p2 != @set p1.z += 0.1
        @test p2 != @set p1.u *= -1
        @test p2 != @set p1.v *= -1
        @test p2 != @set p1.w *= -1
        @test p2 != @set p1.weight += 0.1
        @test p2 != @set p1.new_history = !p1.new_history
    end
end

@testset "validate_record_constants" begin
    R = typeof((E = Float32(1), weight = Float32(2)))
    @test PHSP.validate_record_constants(Bool, R)
    R = typeof((E = Float64(1), weight = Float32(2)))
    @test !PHSP.validate_record_constants(Bool, R)
    R = typeof((E = Float32(1), weight = Float32(2), w=Float32(1)))
    @test !PHSP.validate_record_constants(Bool, R)
end

h = StaticIAEAHeader((E=Float32(1), weight=Float32(2)),0,0)()
@test sizeof(h) === 0

end#module


text:
mutable struct TypeMeta
    group::String
    version:::String
    kind::String
end

mutable struct ObjectMeta
    name::String
    namesapce:::String
    uid::String
    resourceversion::String
    labels::Dict{String, String}
    annotations::Dict{String, String}
end

max_stars_repo_name: pbarker/kubemeta

max_issues_repo_name: pbarker/kubemeta

max_forks_repo_name: pbarker/kubemeta


text:
## Day 18: Snailfish ##
#######################

# Each number is a tree

function reduce!(tree)
    # Search for the leftmost pair to explode
    explode!(tree, 0)
end

# Every number has an index. I think using this index can make the code cleaner.

bigness(tree::Int) = 1
bigness(tree::Array) = sum(bigness.(tree))

@assert bigness([1, 2]) == 2
@assert bigness([[1, 2], 3]) == 3
@assert bigness([9, [8, 7]]) == 3
@assert bigness([[1, 9], [8, 5]]) == 4
@assert bigness([[[[1, 2], [3, 4]], [[5, 6], [7, 8]]], 9]) == 9

function info_about_index(tree::Int, index)
    @assert index == 1
    (tree, 0)
end

function info_about_index(tree::Array, index)
    # Return the number at the given index + information
    left_bigness = 0 # always left_bigness < index
    for t in tree
        sub_bigness = bigness(t)
        if left_bigness + sub_bigness < index
            left_bigness += sub_bigness
        else
            sub_v, sub_d = info_about_index(t, index - left_bigness)
            return (sub_v, sub_d + 1)
        end
    end
    throw("Out of bounds error")
end

@assert info_about_index([1, 2], 1) == (1, 1)
@assert info_about_index([9, [8, 7]], 2) == (8, 2)

function delete_at_index!(tree::Array, index)
    left_bigness = 0
    for i = 1:length(tree)
        t = tree[i]
        sub_bigness = bigness(t)
        if left_bigness + sub_bigness < index
            left_bigness += sub_bigness
        elseif typeof(t) == Int
            popat!(tree, i)
            return tree
        else
            tree[i] = n(t)
            delete_at_index!(tree[i], index - left_bigness)
            if length(tree[i]) == 1
                tree[i] = tree[i][1]
            end
            return tree
        end
    end
    throw("Out of bounds error")
end

@assert [9, 8] == delete_at_index!([9, [8, 7]], 3)
@assert [9, [8, 5], 1] == delete_at_index!([[1, 9], [8, 5], 1], 1)
@assert [1, [8, 5], 1] == delete_at_index!([[1, 9], [8, 5], 1], 2)
@assert [[1, 9], 5, 1] == delete_at_index!([[1, 9], [8, 5], 1], 3)
@assert [[1, 9], 8, 1] == delete_at_index!([[1, 9], [8, 5], 1], 4)
@assert [21, [[0, [7, 6]], [[7, 6], [4, 7]]]] == delete_at_index!(n([21, [[[0, 2], [7, 6]], [[7, 6], [4, 7]]]]), 3)

function set_at!(tree::Array, index, value)
    left_bigness = 0
    for i = 1:length(tree)
        t = tree[i]
        sub_bigness = bigness(t)
        if left_bigness + sub_bigness < index
            left_bigness += sub_bigness
        elseif typeof(t) == Int
            @assert index - left_bigness == 1
            tree[i] = value
            return tree
        else
            set_at!(t, index - left_bigness, value)
            if length(t) == 1
                tree[i] = t[1]
            end
            return tree
        end
    end
    throw("Out of bounds error")
end

@assert [[7, 9], [8, 5], 1] == set_at!([[1, 9], [8, 5], 1], 1, 7)
@assert [[1, 7], [8, 5], 1] == set_at!([[1, 9], [8, 5], 1], 2, 7)
@assert [[1, 9], [7, 5], 1] == set_at!([[1, 9], [8, 5], 1], 3, 7)
@assert [[1, 9], [8, 7], 1] == set_at!([[1, 9], [8, 5], 1], 4, 7)
@assert [[1, 9], [8, 5], 7] == set_at!([[1, 9], [8, 5], 1], 5, 7)

function explode_one!(tree::Array{Any,1})
    for i = 1:bigness(tree)
        left_value, depth = info_about_index(tree, i)
        if depth > 5
            throw("Depth too big, malformed tree")
        elseif depth == 5
            right_value, depth = info_about_index(tree, i + 1)

            if i - 1 >= 1
                set_at!(tree, i - 1, left_value + info_about_index(tree, i - 1)[1])
            end
            if i + 2 <= bigness(tree)
                set_at!(tree, i + 2, right_value + info_about_index(tree, i + 2)[1])
            end
            set_at!(tree, i, 0)
            delete_at_index!(tree, i + 1)
            return tree
        end
    end
end

"""Normalize away the type issuses I see when entering literals like [[1,2],[3,4]]"""
function n(tree::Array)::Array{Any,1}
    result = Array{Any,1}()
    for t in tree
        push!(result, n(t))
    end
    tree
end

n(tree::Int)::Int = tree

@assert [1, [8, 5]] == delete_at_index!(n([[1, 9], [8, 5]]), 2)

@assert [[[[0, 9], 2], 3], 4] == explode_one!([[[[[9, 8], 1], 2], 3], 4])
@assert [7, [6, [5, [7, 0]]]] == explode_one!([7, [6, [5, [4, [3, 2]]]]])
@assert [[6, [5, [7, 0]]], 3] == explode_one!([[6, [5, [4, [3, 2]]]], 1])
@assert [[3, [2, [8, 0]]], [9, [5, [4, [3, 2]]]]] == explode_one!(n([[3, [2, [1, [7, 3]]]], [6, [5, [4, [3, 2]]]]]))
@assert [[3, [2, [8, 0]]], [9, [5, [7, 0]]]] == explode_one!(n([[3, [2, [8, 0]]], [9, [5, [4, [3, 2]]]]]))

function split_one!(tree::Array{Any,1})
    for i = 1:length(tree)
        t = tree[i]
        if typeof(t) == Int && t >= 10
            tree[i] = n([Int(floor(t / 2)), Int(ceil(t / 2))])
            return (tree, :split)
        elseif typeof(t) == Int
            continue
        else
            tree[i] = n(t)
            _, action = split_one!(tree[i])
            if action == :split
                return (tree, :split)
            end
        end
    end
    (tree, :no_split)
end

@assert ([[[[0, 7], 4], [[7, 8], [0, 13]]], [1, 1]], :split) == split_one!(n([[[[0, 7], 4], [15, [0, 13]]], [1, 1]]))
@assert ([[[[0, 7], 4], [[7, 8], [0, [6, 7]]]], [1, 1]], :split) == split_one!(n([[[[0, 7], 4], [[7, 8], [0, 13]]], [1, 1]]))
@assert ([[[[0, 7], 4], [[7, 8], [0, [6, 7]]]], [1, 1]], :no_split) == split_one!(n([[[[0, 7], 4], [[7, 8], [0, [6, 7]]]], [1, 1]]))

function reduce!(tree::Array{Any,1})
    while true
        # Try to explode a pair
        old_bigness = bigness(tree)
        explode_one!(tree)
        if bigness(tree) != old_bigness
            continue
        end

        # Try to split a pair
        _, action = split_one!(tree)
        if action == :split
            continue
        end

        # No more changes
        return tree
    end
end

@assert [[[[0, 7], 4], [[7, 8], [6, 0]]], [8, 1]] == reduce!(n([[[[[4, 3], 4], 4], [7, [[8, 4], 9]]], [1, 1]]))

include("test-17.jl")
include("input-17.jl")

function final_sum(list_of_trees)
    sum = list_of_trees[1]
    for i = 2:length(list_of_trees)
        sum = reduce!(n([sum, deepcopy(list_of_trees[i])]))
    end
    sum
end

@assert [[[[6, 6], [7, 6]], [[7, 7], [7, 0]]], [[[7, 7], [7, 7]], [[7, 8], [9, 9]]]] == final_sum(day17_test)

magnitude(tree::Int) = tree
magnitude(tree::Array) = 3 * magnitude(tree[1]) + 2 * magnitude(tree[2])

@assert 4140 == magnitude([[[[6, 6], [7, 6]], [[7, 7], [7, 0]]], [[[7, 7], [7, 7]], [[7, 8], [9, 9]]]])
@assert 4469 == magnitude(final_sum(day17_input))

include("test-17.jl")
include("input-17.jl")

function max_magnitude(list_of_trees)
    best = 0
    for i = 1:length(list_of_trees), j = 1:length(list_of_trees)
        if i == j
            continue
        end
        best = max(best, magnitude(final_sum(n(deepcopy([list_of_trees[i], list_of_trees[j]])))))
    end
    best
end

@assert 3993 == max_magnitude(day17_test)
@assert 4770 == max_magnitude(day17_input)


text:
import Random

perturbation(step, sig) = -sig + 2*sig*rand()

function initial_point(N, side, tol; seed=123)

	if seed > 0
		Random.seed!(seed)
	else
		Random.seed!()
	end

	p = [ zeros(2) for i in 1:N ]

	l = ceil(Int, sqrt(N))
	step = (side-tol)/l

	if tol > step
		error("tol must be smaller than $step")
	end

	skip = Int[]

	while length(skip) < (l^2-N)

		i = rand(1:l^2)

		if findfirst(isequal(i), skip) == nothing
			push!(skip, i)
		end
	end

	sig = (step-tol)/2

	ip = 0
	igrid = 0
	x = [ tol/2, tol/2 ]

	for i in 1:l
		for j in 1:l

			igrid += 1

			if findfirst(isequal(igrid), skip) != nothing
				continue
			end

			ip += 1
			p[ip] .= x .+ perturbation(step, sig)
			x[2] += step
		end

		x[2] = tol/2
		x[1] += step
	end

	return p
end

initial_point(data :: Data; tol=0.9, seed=123) = initial_point(data.N, data.side, tol, seed=seed)


text:
using NCDatasets

project_root_dir = pwd()
data_dir = joinpath(project_root_dir)
domain_dir = joinpath(project_root_dir)
casename = "Sandbox"

forcing_file = joinpath(data_dir, "ocn_forcing.nc")

config = Dict{Any, Any}(

    "DRIVER" => Dict(
        "casename"           => casename,
        "caseroot"           => joinpath(project_root_dir, casename, "caseroot"),
        "caserun"            => joinpath(project_root_dir, casename, "caserun"),
        "archive_root"       => joinpath(project_root_dir, casename, "archive"),
    ),

    "MODEL_MISC" => Dict(
        "timetype"               => "DateTimeNoLeap",
        "init_file"              => "",
        "rpointer_file"          => "rpointer.iom",
        "daily_record"           => [],
        "monthly_record"         => ["{ESSENTIAL}",],
        "enable_archive"         => true,
    ),

    "MODEL_CORE" => Dict(

        "domain_file"                  => joinpath(domain_dir, "domain.nc"),
        "topo_file"                    => "",

        "cdata_var_file_map"           => Dict(
            "HMXL"      => forcing_file,
            "TEMP"      => forcing_file,
            "SALT"      => forcing_file,
            "QFLX_TEMP" => forcing_file,
            "QFLX_SALT" => forcing_file,
        ),

        "cdata_beg_time"               => "0001-01-01 00:00:00",
        "cdata_end_time"               => "0002-01-01 00:00:00",
        "cdata_align_time"             => "0001-01-01 00:00:00",

        "z_w"                          => nothing,

        "substeps"                     => 8,
        "MLD_scheme"                   => "datastream",
        "Qflx"                         => "on",
        "Qflx_finding"                 => "off",
        "convective_adjustment"        => "on",
        "advection_scheme"             => "ekman_AGA2020_allowU",

        "weak_restoring"               => "on",
        "τwk_TEMP"                     => 86400.0 * 365 * 1000,
        "τwk_SALT"                     => 86400.0 * 365 * 1000,

        "τ_frz"                        => 3600.0,
        "Ekman_layers"                 => 2,
        "Returnflow_layers"            => 8,

        "transform_vector_field"       => true,
    ),

)

Dataset("ocn_forcing.nc", "r") do ds

    z_w_top = nomissing(ds["z_w_top"][:], NaN)
    z_w_bot = nomissing(ds["z_w_bot"][:], NaN)

    z_w      = zeros(Float64, length(z_w_top)+1)
    z_w[1:end-1]   = z_w_top
    z_w[end] = z_w_bot[end]

    global config["MODEL_CORE"]["z_w"]  = z_w
end

using TOML

open("config.toml", "w") do io
    TOML.print(io, config)
end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# # Cohomology, Homology, and Representatives

# In this section, we will show how Ripserer can be used to find critical simplices and
# representative (co)cycles.

# We start by loading some packages and generating some data.

using LinearAlgebra
using Plots
using Ripserer
using Random # hide
Random.seed!(1337) # hide
gr() # hide
nothing # hide

function annulus(n, r1=1, r2=2, offset=(0, 0))
    result = Tuple{Float64,Float64}[]
    while length(result) < n
        point = 2 * r2 * rand(2) .- r2
        if r1 < norm(point) < r2
            push!(result, (point[1] + offset[1], point[2] + offset[2]))
        end
    end
    return result
end

data = annulus(300)

scatter(data; label="data", markersize=2, aspect_ratio=1)

# Let's start by taking a look at the persistence diagram.

diagram = ripserer(data)
plot(diagram)

# The diagram tells us that there is a persistent hole in the data, but tells us nothing
# about the location of the hole. Ripserer provides several methods to locate it. We'll
# start with the simplest.

# ## Critical simplices

# The first option is to find the death simplex of the interval.

# We start by extracting the interval in question. Keep in mind that the diagrams are sorted
# by persistence, so the last element will always be the most persistent.

most_persistent = diagram[2][end]

# Notice that the interval has two simplices attached to it, the birth simplex and the death
# simplex. We can extract them with `birth_simplex` and `death_simplex` respectively.

death_sx = death_simplex(most_persistent)

# A simplex acts just like an array of indices, so it can be used to index into the data.

data[death_sx]

# Ripserer also provides a Plots recipe for plotting simplices. It is invoked by passing the
# simplex and the data to `plot`. Not that only the edges of the simplices are plotted.

scatter(data; label="data", markersize=2, aspect_ratio=1)
plot!(death_sx, data; label="death simplex")
plot!(birth_simplex(diagram[2][end]), data; label="birth simplex")

# The birth simplex is the simplex that first connects the hole. The death simplex is the
# simplex that fills the hole in.

# While the death simplex gives us a vague idea of where the hole is located, there are
# other methods available.

# ## Representative Cocycles

# By default, Ripserer computes persistent cohomology. The resulting diagrams of persistent
# homology and cohomology are the same, but computing cohomology is much more
# efficient. When computing persistent cohomology, we can tell Ripserer to also compute
# representative cocycles. This is controlled with the `reps` keyword argument.

# Let's take a look at the most persistent cocycle of our data set.

diagram_cocycles = ripserer(data; reps=true)
most_persistent_co = diagram_cocycles[2][end]

# Notice that now, the interval also has a `representative` attached. The representative is
# an array of pairs `Simplex => value`, where the value is the coefficient of the simplex.
# In reality, the type is different, but it acts exactly the same as a `Pair`.

cocycle = representative(most_persistent_co)

# The representative can be plotted in the same way as a simplex.

scatter(data; label="data", markersize=2, aspect_ratio=1)
plot!(cocycle, data; label="cocycle")

# The cocycle is a collection of 1-simplices that, if removed, would break the cycle in our
# data set. This does not correspond to most people's intuitive understanding of a hole, but
# it can be useful in some situations. To find something more intuitive, we have to look to
# homology and its representative cycles.

# ## Representative Cycles

# Ripserer supports two algorithms for computing representative cocycles. One is computing
# persistent homology directly, and the other is involuted homology computation. Involuted
# homology computes cohomology first and then uses its result to recompute cycles. While
# this increases the running time somewhat, it is still usually much more efficient than
# computing persistent homology directly. The difference is especially large for filtrations
# where the number of simplices increases quickly with dimension, such as Vietoris-Rips
# filtrations. See [this paper](https://arxiv.org/abs/2105.03629) for more information.

# Involuted homology is computed by passing the argument `alg=:involuted` to `ripserer`. If
# we wanted direct homology computation, we would use `alg=:homology`. The results for both
# cases are exactly the same.

# Let's try it out. Note that invoking homology also turns on `reps` for dimensions one and
# higher.

diagram_cycles = ripserer(data; alg=:involuted)
most_persistent_ho = diagram_cycles[2][end]

# If an interval with a representative is passed to `plot`, the representative is plotted.

scatter(data; label="data", markersize=2, aspect_ratio=1)
plot!(most_persistent_ho, data; label="cycle")

# The cycle is still not the prettiest, but it at least corresponds to a topological circle
# wound around the hole in the middle of the data set. Sometimes, the cycle will also have
# multiple connected components. All except one will be contractible at the time the cycle
# exists. To make the result look even better, we can try reconstructing the shortest
# representative cycle.

# ## Reconstructed Shortest Cycles

# This method uses information from the cocycle to reconstruct the shortest
# cycle. Essentially what the method does is that it picks an edge in the cocycle and
# connects it through the edges in the filtration at the specified time. A limitation of
# this method is that it only works for one-dimensional representatives.

# Reconstruction is done on a per interval basis, as it would take a long time to
# reconstruct all cycles in large diagrams.

# Let's start with a basic reconstruction.

filtration = diagram_cocycles[2].filtration
reconstructed_at_birth = reconstruct_cycle(filtration, most_persistent_co)

scatter(data; label="data", markersize=2, aspect_ratio=1)
plot!(reconstructed_at_birth, data; label="reconstruction")

# This looks much nicer than the homology example, but could still use some improvement. To
# improve it, we can set a time at which to reconstruct the cycle.

# As time goes on and more simplices are added to the filtration, the shapes of the shortest
# cycles change as well. The previous example was drawn at interval birth time, which is the
# default. Let's see what happens if we set the time to the interval midpoint.

midpoint = (death(most_persistent_co) - birth(most_persistent_co)) / 2
reconstructed_at_midpoint = reconstruct_cycle(filtration, most_persistent_co, midpoint)

scatter(data; label="data", markersize=2, aspect_ratio=1)
plot!(reconstructed_at_midpoint, data; label="reconstruction")

# As an extreme case, let's look at what the cycle looks like right before its death.

scatter(data; label="data", markersize=2, aspect_ratio=1)
plot!(
    reconstruct_cycle(filtration, most_persistent_co, death(most_persistent_co) - 0.01),
    data;
    label="reconstruction",
)

# As the time nears the death time, the cycle gets closer to looking like the death simplex.


text:
using OPCSPs.MVNTools

using PyCall
using PyPlot
@pyimport matplotlib.patches as patches

function plot_confidence_ellipse(d::MVN, num_sigma::Int=3)
    assert(length(d.mean)==2)
    lambda, phi = eig(d.covariance)
    ell = patches.Ellipse(xy = d.mean,
                          width = sqrt(lambda[1])*num_sigma*2,
                          height = sqrt(lambda[2])*num_sigma*2,
                          angle = 180/pi*atan2(phi[2,1], phi[1,1]))
    ell[:set_facecolor]("none")
    gca()[:add_artist](ell)
    draw()
end

d = MVN([0,1], [1 0.08; 0.08 0.04])

n_samples = 10000
samples = Array(Float64, n_samples, 2)
rng = MersenneTwister()
for i in 1:n_samples
    samples[i,:] = rand!(rng, Array(Float64, length(d)), d)
end

figure(1)
clf()
plot_confidence_ellipse(d, 1)
plot_confidence_ellipse(d, 2)
plot_confidence_ellipse(d, 3)
xlim(-3, 3)
ylim(-3, 3)
scatter(samples[:,1], samples[:,2])


text:
using Documenter, Bigsimr
using Distributions

DocMeta.setdocmeta!(Bigsimr, :DocTestSetup, :(using Bigsimr); recursive=true)

makedocs(
    sitename = "Bigsimr.jl",
    modules  = [Bigsimr],
    doctest  = false,
    pages    = [
        "Bigsimr.jl" => "index.md",
        "Guides" => [
            "Getting Started" => "getting_started.md",
            "Pearson Matching" => "pearson_matching.md",
            "Nearest Correlation Matrix" => "nearest_correlation_matrix.md"
        ],
        "API Reference" => "main_functions.md",
        "Details" => "details.md",
        "Index" => "function_index.md"
    ]
)

deploydocs(
    repo = "github.com/SchisslerGroup/Bigsimr.jl.git",
    versions = ["stable" => "v^", "v#.#", "dev" => "master"]
)