text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami<nakagami@gmail.com>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################

using Sockets

const PLUGIN_LIST = "Srp256,Srp"
const BUFFER_LEN = 1024
const MAX_CHAR_LENGTH = 32767
const BLOB_SEGMENT_SIZE = 32000

function DEBUG_OUTPUT(s)
    # println(s)
end

function INFO_SQL_SELECT_DESCRIBE_VARS()::Vector{UInt8}
    [
        isc_info_sql_select,
        isc_info_sql_describe_vars,
        isc_info_sql_sqlda_seq,
        isc_info_sql_type,
        isc_info_sql_sub_type,
        isc_info_sql_scale,
        isc_info_sql_length,
        isc_info_sql_null_ind,
        isc_info_sql_field,
        isc_info_sql_relation,
        isc_info_sql_owner,
        isc_info_sql_alias,
        isc_info_sql_describe_end
    ]
end

mutable struct WireChannel
    socket::Union{TCPSocket, Nothing}
    in::Union{Arc4, ChaCha20, Nothing}
    out::Union{Arc4, ChaCha20, Nothing}
    function WireChannel(host::String, port::UInt16)
        socket = Sockets.connect(host, port)
        new(socket, nothing, nothing)
    end
end

function set_crypt_key(chan::WireChannel, key::Vector{UInt8}, plugin::String, nonce::Vector{UInt8})
    if plugin == "Arc4"
        chan.in = Arc4(key)
        chan.out = Arc4(key)
    elseif plugin == "ChaCha"
        key = sha256(key)
        chan.in = ChaCha20(key, nonce, Int32(0))
        chan.out = ChaCha20(key, nonce, Int32(0))
    end
end

function recv(chan::WireChannel, nbytes::Int)::Vector{UInt8}
    data::Vector{UInt8} = zeros(UInt8, nbytes)
    read!(chan.socket, data)
    if chan.in != nothing
        data = translate(chan.in, data)
    end
    data
end

function send(chan::WireChannel, data::Vector{UInt8})
    if chan.out != nothing
        data = translate(chan.out, data)
    end
    write(chan.socket, data)
end

function close!(chan::WireChannel)
    if chan.socket != nothing
        close(chan.socket)
        chan.socket = nothing
    end
end

function isopen(chan::WireChannel)::Bool
    chan.socket != nothing
end

mutable struct WireProtocol
    buf::Vector{UInt8}

    chan::WireChannel
    host::String
    port::UInt16
    username::String
    password::String

    db_handle::Int32
    trans_handle::Int32

    protocol_version::Int32
    accept_architecture::Int32
    accept_type::Int32
    lazy_response_count::Int

    accept_plugin_name::String
    auth_data::Vector{UInt8}

    timezone::String

    function WireProtocol(host::AbstractString, user::AbstractString, password::AbstractString, port::UInt16)
        chan = WireChannel(host, port)
        new([], chan, host, port, user, password, -1, -1, -1, -1, -1, 0, "", [], "")
    end
end

function pack_uint32(wp::WireProtocol, i::Int)
    wp.buf = vcat(wp.buf, UInt8[UInt8(i >> 24 & 0xFF), UInt8(i >> 16 & 0xFF), UInt8(i >> 8 & 0xFF), UInt8(i & 0xFF)])
end

function pack_uint32(wp::WireProtocol, i::Int32)
    wp.buf = vcat(wp.buf, UInt8[UInt8(i >> 24 & 0xFF), UInt8(i >> 16 & 0xFF), UInt8(i >> 8 & 0xFF), UInt8(i & 0xFF)])
end

function pack_uint32(wp::WireProtocol, i::UInt32)
    # pack big endian uint32
    wp.buf = vcat(wp.buf, UInt8[UInt8(i >> 24 & 0xFF), UInt8(i >> 16 & 0xFF), UInt8(i >> 8 & 0xFF), UInt8(i & 0xFF)])
end

function pack_bytes(wp::WireProtocol, b::Vector{UInt8})::Vector{UInt8}
    wp.buf = vcat(wp.buf, xdr_bytes(b))
end

function pack_string(wp::WireProtocol, s::String)
    wp.buf = vcat(wp.buf, xdr_bytes(Vector{UInt8}(s)))
end

function append_bytes(wp::WireProtocol, b::Vector{UInt8})
    wp.buf = vcat(wp.buf, b)
end

function get_srp_client_public_bytes(client_public::BigInt)::Vector{UInt8}
    b::Vector{UInt8} = Vector{UInt8}(bytes2hex(bigint_to_bytes(client_public)))
    if length(b) > 254
        vcat(
            UInt8[CNCT_specific_data, 255, 0],
            b[1:254],
            UInt8[CNCT_specific_data, length(b) - 254 + 1, 1],
            b[255:length(b)]
        )
    else
        vcat(
            UInt8[CNCT_specific_data, length(b) + 1, 0],
            b[1:length(b) + 1],
            UInt8[CNCT_specific_data, length(b) - 254 + 1, 1],
            b
        )
    end
end

function uid(username::String, password::String, auth_plugin_name::String, wire_crypt::Bool, client_public::BigInt)::Vector{UInt8}
    sys_user = if haskey(ENV, "USER")
            ENV["USER"]
        elseif haskey(ENV, "USERNAME")
            ENV["USERNAME"]
        else
            ""
        end
    sys_user_bytes = Vector{UInt8}(sys_user)
    hostname_bytes = Vector{UInt8}(gethostname())
    plugin_list_name_bytes = Vector{UInt8}(PLUGIN_LIST)
    plugin_name_bytes = Vector{UInt8}(auth_plugin_name)
    username_bytes = Vector{UInt8}(username)
    specific_data = get_srp_client_public_bytes(client_public)
    vcat(
        UInt8[CNCT_login, length(username_bytes)], username_bytes,
        UInt8[CNCT_plugin_name, length(plugin_name_bytes)], plugin_name_bytes,
        UInt8[CNCT_plugin_list, length(plugin_list_name_bytes)], plugin_list_name_bytes,
        specific_data,
        UInt8[CNCT_client_crypt, 4, wire_crypt ? 1 : 0, 0, 0, 0],
        UInt8[CNCT_user, length(sys_user_bytes)], sys_user_bytes,
        UInt8[CNCT_host, length(hostname_bytes)], hostname_bytes,
        UInt8[CNCT_user_verification, 0]
    )
end

function send_packets(wp::WireProtocol)
    send(wp.chan, wp.buf)
    wp.buf = []
end

function suspend_buffer(wp::WireProtocol)::Vector{UInt8}
    deepcopy(wp.buf)
end

function resume_buffer(wp::WireProtocol, buf::Vector{UInt8})
    append!(wp.buf, buf)
end

function recv_packets(wp::WireProtocol, n::Int)::Vector{UInt8}
    recv(wp.chan, n)
end

function recv_packets(wp::WireProtocol, n::UInt32)::Vector{UInt8}
    recv_packets(wp, Int(n))
end

function recv_packets_alignment(wp::WireProtocol, n::UInt32)::Vector{UInt8}
    padding = n % 4
    if padding > 0
        padding = 4 - padding
    end
    buf = recv_packets(wp, UInt32(n + padding))
    buf[1:n]
end

function recv_packets_alignment(wp::WireProtocol, n::Int32)::Vector{UInt8}
    recv_packets_alignment(wp, UInt32(n))
end

function recv_packets_alignment(wp::WireProtocol, n::Int64)::Vector{UInt8}
    recv_packets_alignment(wp, UInt32(n))
end

function parse_status_vector(wp::WireProtocol)::Tuple{Vector{UInt32}, Int, String}
    sql_code::Int = 0
    gds_code = 0
    gds_codes::Vector{UInt32} = []
    num_arg = 0
    message::String = ""

    errmsgs = get_errmsgs()

    n = bytes_to_buint32(recv_packets(wp, 4))
    while n != isc_arg_end
        if n == isc_arg_gds
            gds_code = bytes_to_buint32(recv_packets(wp, 4))
            if gds_code != 0
                push!(gds_codes, gds_code)
                message *= errmsgs[gds_code]
                num_arg = 0
            end
        elseif n == isc_arg_number
            num = bytes_to_buint32(recv_packets(wp, 4))
            if gds_code == 335544436
                sql_code = num
            end
            num_arg += 1
            message = replace.(message, [string("@", num_arg)=>num])[1]
        elseif n == isc_arg_string
            nbytes = bytes_to_buint32(recv_packets(wp, 4))
            s = String(recv_packets_alignment(wp, nbytes))
            num_arg += 1
            message = replace.(message, [string("@", num_arg)=>s])[1]
        elseif n == isc_arg_iterpreted
            nbytes = bytes_to_buint32(recv_packets(wp, 4))
            message *= String(recv_packets_alignment(wp, nbytes))
        elseif n == isc_arg_sql_state
            nbytes = bytes_to_buint32(recv_packets(wp, 4))
            recv_packets_alignment(wp, nbytes)  # skip status code
        end
        n = bytes_to_buint32(recv_packets(wp, 4))
    end

    (gds_codes, sql_code, message)
end

function parse_op_response(wp::WireProtocol)::Tuple{Int32, Vector{UInt8}, Vector{UInt8}}
    h = bytes_to_bint32(recv_packets(wp, 4))            # Object handle
    oid = recv_packets(wp, 8)                           # Object ID
    buf_len = Int(bytes_to_bint32(recv_packets(wp, 4))) # buffer length
    buf = recv_packets_alignment(wp, buf_len)

    gds_codes, sql_code, message = parse_status_vector(wp)
    if length(gds_codes) != 0
        throw(DomainError("response error", message))
    end
    (h, oid, buf)
end

function guess_wire_crypt(buf::Vector{UInt8})::Tuple{String, Vector{UInt8}}
    params = Dict()
    i = 1
    while i <= length(buf)
        k = buf[i]
        i += 1
        ln = buf[i]
        i += 1
        v = buf[i:i+ln-1]
        i += ln
        params[k] = v
    end
    if haskey(params, 3) && params[3][1:7] == b"ChaCha\x00"
        return "ChaCha", params[3][8:length(params[3])-4]
    end

    "Arc4", []
end

function parse_connect_response(wp::WireProtocol, username::String, password::String, wire_crypt::Bool, client_public::BigInt, client_secret::BigInt)
    DEBUG_OUTPUT("parse_connect_response")
    op_code = bytes_to_bint32(recv_packets(wp, 4))
    while op_code == op_dummy
        op_code = bytes_to_bint32(recv_packets(wp, 4))
    end
    if op_code == op_reject
        throw(DomainError("op_reject"))
    end
    if op_code == op_response
        parse_op_response(wp)
        # NOT REACH HERE?
        throw(DomainError("op_response"))
    end

    wp.protocol_version = recv_packets(wp, 4)[4]
    wp.accept_architecture = bytes_to_bint32(recv_packets(wp, 4))
    wp.accept_type = bytes_to_bint32(recv_packets(wp, 4))

    @assert op_code == op_cond_accept || op_code == op_accept_data

    ln = bytes_to_buint32(recv_packets(wp, 4))
    data = recv_packets_alignment(wp, ln)

    ln = bytes_to_buint32(recv_packets(wp, 4))
    wp.accept_plugin_name = String(recv_packets_alignment(wp, ln))

    # is_authenticated == 0
    @assert bytes_to_buint32(recv_packets(wp, 4)) == 0

    # skip keys
    ln = bytes_to_buint32(recv_packets(wp, 4))
    recv_packets_alignment(wp, ln)

    @assert wp.accept_plugin_name == "Srp" || wp.accept_plugin_name == "Srp256"

    if length(data) == 0
        _op_cont_auth(wp, bigint_to_bytes(client_public), wp.accept_plugin_name, "")
        @assert bytes_to_bint32(recv_packets(wp, 4)) == op_cont_auth

        ln = bytes_to_buint32(recv_packets(wp, 4))
        data = recv_packets_alignment(wp, ln)

        # skip plugin name
        ln = bytes_to_buint32(recv_packets(wp, 4))
        recv_packets_alignment(wp, ln)

        # skip plugin name list
        ln = bytes_to_buint32(recv_packets(wp, 4))
        recv_packets_alignment(wp, ln)

        # skip keys
        ln = bytes_to_buint32(recv_packets(wp, 4))
        recv_packets_alignment(wp, ln)
    end
    ln = bytes_to_uint16(data[1:2])
    server_salt = data[3:ln+2]
    server_public_string = data[5+ln:length(data)]
    if length(server_public_string) % 2 != 0
        server_public_string = vcat(Vector{UInt8}([0x30]), server_public_string)
    end
    server_public = bytes_to_bigint(hex2bytes(server_public_string))

    auth_data, session_key = get_client_proof(
        uppercase(username),
        password,
        server_salt,
        client_public,
        server_public,
        client_secret,
        wp.accept_plugin_name,
    )
    if op_code == op_cond_accept
        _op_cont_auth(wp, auth_data, wp.accept_plugin_name, "")
        _, _, buf = _op_response(wp)
        encrypt_plugin, nonce = guess_wire_crypt(buf)
    end

    if wire_crypt && session_key != nothing
        _op_crypt(wp, encrypt_plugin)
        set_crypt_key(wp.chan, session_key, encrypt_plugin, nonce)
        _op_response(wp)
    else
        wp.auth_data = auth_data
    end

end

function parse_select_items(wp::WireProtocol, buf::Vector{UInt8}, xsqlda::Vector{XSQLVAR})::Int
    ln = 0
    index = 0
    i = 1

    item = buf[i]
    while i <= length(buf) && buf[i] != isc_info_end
        item = buf[i]
        i += 1
        if item ==  isc_info_sql_sqlda_seq
            ln = bytes_to_int16(buf[i:i+1])
            i += 2
            index = bytes_to_int32(buf[i:i+ln-1])
            i += ln
        elseif item == isc_info_sql_type
            ln = bytes_to_int16(buf[i:i+1])
            i += 2
            sqltype = bytes_to_int32(buf[i:i+ln-1])
            if sqltype % 2 != 0
                    sqltype -= 1
            end
            xsqlda[index].sqltype = sqltype
            i += ln
        elseif item == isc_info_sql_sub_type
            ln = bytes_to_int16(buf[i:i+1])
            i += 2
            xsqlda[index].sqlsubtype = bytes_to_int32(buf[i:i+ln-1])
            i += ln
        elseif item == isc_info_sql_scale
            ln = bytes_to_int16(buf[i:i+1])
            i += 2
            xsqlda[index].sqlscale = bytes_to_int32(buf[i:i+ln-1])
            i += ln
        elseif item == isc_info_sql_length
            ln = bytes_to_int16(buf[i:i+1])
            i += 2
            xsqlda[index].sqllen = bytes_to_int32(buf[i:i+ln-1])
            i += ln
        elseif item == isc_info_sql_null_ind
            ln = bytes_to_int16(buf[i:i+1])
            i += 2
            xsqlda[index].null_ok = bytes_to_int32(buf[i:i+ln-1]) != 0
            i += ln
        elseif item == isc_info_sql_field
            ln = bytes_to_int16(buf[i:i+1])
            i += 2
            xsqlda[index].fieldname = bytes_to_str(buf[i:i+ln-1])
            i += ln
        elseif item == isc_info_sql_relation
            ln = bytes_to_int16(buf[i:i+1])
            i += 2
            xsqlda[index].relname = bytes_to_str(buf[i:i+ln-1])
            i += ln
        elseif item == isc_info_sql_owner
            ln = bytes_to_int16(buf[i : i+1])
            i += 2
            xsqlda[index].ownname = bytes_to_str(buf[i:i+ln-1])
            i += ln
        elseif item == isc_info_sql_alias
            ln = bytes_to_int16(buf[i : i+1])
            i += 2
            xsqlda[index].aliasname = bytes_to_str(buf[i:i+ln-1])
            i += ln
        elseif item == isc_info_truncated
            return index        # return next index
        elseif item == isc_info_sql_describe_end
            # NOTHING
        else
            throw(DomainError("Invalid item", "$i:$(buf[i])"))
        end
    end

    -1  # no more info
end

function parse_xsqlda(wp::WireProtocol, buf::Vector{UInt8}, stmt_handle::Int32)::Tuple{Int, Vector{XSQLVAR}}
    stmt_type::Int = 0
    xsqlda::Vector{XSQLVAR} = []
    i = 1

    while i <= length(buf)
        if buf[i] == UInt8(isc_info_sql_stmt_type) && buf[i+1] == UInt8(0x04) && buf[i+2] == UInt8(0x00)
            i += 1
            ln = bytes_to_int16(buf[i : i+1])
            i += 2
            stmt_type = bytes_to_int32(buf[i : i+ln-1])
            i += ln
        elseif buf[i] == UInt8(isc_info_sql_select) && buf[i+1] == UInt8(isc_info_sql_describe_vars)
            i += 2
            ln = bytes_to_int16(buf[i : i+1])
            i += 2
            col_len = bytes_to_int32(buf[i:i+ln-1])
            if col_len != 0
                for _ in range(1, length=col_len)
                    xsqlvar = XSQLVAR(0, 0, 0, 0, false, "", "", "", "")
                    push!(xsqlda, xsqlvar)
                end
                next_index = parse_select_items(wp, buf[i+ln:length(buf)-1], xsqlda)
                while next_index > 0    # more describe vars
                    _op_info_sql(stmt_handle,
                        vcat(
                            Vector{UInt8}([isc_info_sql_sqlda_start, 2]),
                            int16_to_bytes(int16(next_index)),
                            INFO_SQL_SELECT_DESCRIBE_VARS(),
                        )
                    )
                    _, _, buf = op_response(wp)
                    # buf[1:2] == [0x04,0x07]
                    ln = bytes_to_int16(buf[3:4])
                    # bytes_to_int(buf[5:5+ln]) == col_len
                    next_index = p._parse_select_items(buf[5+ln:length(buf)-1], xsqlda)
                end
            end
        else
            break
        end
    end

    stmt_type, xsqlda
end

function get_blob_segments(wp::WireProtocol, blob_id::Vector{UInt8}, trans_handle::Int32)::Vector{UInt8}
    suspend_buf = suspend_buffer(wp)
    blob::Vector{UInt8} = []
    _op_open_blob(wp, blob_id, trans_handle)
    blob_handle, _, _ = _op_response(wp)

    more_data = 1
    while more_data != 2
        _op_get_segment(wp, blob_handle)
        more_data, _, buf = _op_response(wp)
        while length(buf) > 0
            ln = bytes_to_int16(buf[1:2])
            blob = vcat(blob, buf[3:ln+2])
            buf = buf[ln+3:length(buf)]
        end
    end

    _op_close_blob(wp, blob_handle)
    if wp.accept_type == ptype_lazy_send
        wp.lazy_response_count += 1
    else
        _op_response(wp)
    end

    resume_buffer(wp, suspend_buf)
    blob
end

function _op_connect(wp::WireProtocol, db_name::String, username::String, password::String, wire_crypt::Bool, client_public::BigInt)
    DEBUG_OUTPUT("_op_connect")
    # PROTOCOL_VERSION, Arch type (Generic=1), min, max, weight
    protocols = hex2bytes(string(
        "ffff800d00000001000000000000000500000008", # 13, 1, 0, 5, 8
        "ffff800e0000000100000000000000050000000a", # 14, 1, 0, 5, 10
        "ffff800f0000000100000000000000050000000c", # 15, 1, 0, 5, 12
        "ffff80100000000100000000000000050000000e", # 16, 1, 0, 5, 14
        "ffff801100000001000000000000000500000010", # 17, 1, 0, 5, 16
    ))
    protocols_len = div(length(protocols), 20)

    pack_uint32(wp, op_connect)
    pack_uint32(wp, op_attach)
    pack_uint32(wp, 3)  # CONNECT_VERSION3
    pack_uint32(wp, 1)  # Arc type(GENERIC)
    pack_string(wp, db_name)
    pack_uint32(wp, protocols_len)  # number of protocols
    pack_bytes(wp, uid(username, password, "Srp256", wire_crypt, client_public))
    append_bytes(wp, protocols)
    send_packets(wp)
end

function _op_create(wp::WireProtocol, db_name::String, username::String, password::String, timezone::String, page_size::Int32)
    DEBUG_OUTPUT("_op_create")
    encode = b"UTF8"

    username_bytes = Vector{UInt8}(username)
    password_bytes = Vector{UInt8}(password)
    dpb::Vector{UInt8} = []

    dpb = vcat(
        [isc_dpb_version1],
        [isc_dpb_set_db_charset, UInt8(length(encode))], encode,
        [isc_dpb_lc_ctype, UInt8(length(encode))], encode,
        [isc_dpb_user_name, UInt8(length(username_bytes))], username_bytes,
        [isc_dpb_password, UInt8(length(password_bytes))], password_bytes,
        [isc_dpb_sql_dialect, 4], int32_to_bytes(Int32(3)),
        [isc_dpb_force_write, 4], int32_to_bytes(Int32(1)),
        [isc_dpb_overwrite, 4], int32_to_bytes(Int32(1)),
        [isc_dpb_page_size, 4], int32_to_bytes(page_size),
    )
    if timezone != ""
        timezone_bytes = Vector{UInt8}(timezone)
        dpb = vcat(dpb, [isc_dpb_session_time_zone, UInt8(length(timezone_bytes))], timezone_bytes)
    end

    if length(wp.auth_data) != 0
        specific_auth_data = Vector{UInt8}(bytes2hex(wp.auth_data))
        dpb = vcat(dpb, [isc_dpb_specific_auth_data, length(specific_auth_data)], specific_auth_data)
    end

    if wp.timezone != ""
        tzname_bytes = Vector{UInt8}(wp.timezone)
        dpb = vcat(dpb, [isc_dpb_session_time_zone, UInt8(length(timezone_bytes))], timezone_bytes)
    end

    pack_uint32(wp, op_create)
    pack_uint32(wp, 0)
    pack_string(wp, db_name)
    pack_bytes(wp, dpb)
    send_packets(wp)
end

function _op_attach(wp::WireProtocol, db_name::String, username::String, password::String, timezone::String)
    DEBUG_OUTPUT("_op_attach")
    encode = b"UTF8"

    username_bytes = Vector{UInt8}(username)
    password_bytes = Vector{UInt8}(password)

    dpb::Vector{UInt8} = vcat(
        [UInt8(isc_dpb_version1)],
        [UInt8(isc_dpb_set_db_charset), UInt8(length(encode))], encode,
        [UInt8(isc_dpb_lc_ctype), UInt8(length(encode))], encode,
        [UInt8(isc_dpb_user_name), UInt8(length(username_bytes))], username_bytes,
        [UInt8(isc_dpb_password), UInt8(length(password_bytes))], password_bytes,
        [UInt8(isc_dpb_sql_dialect), UInt8(4)], int32_to_bytes(Int32(3)),
    )

    if timezone != ""
        timezone_bytes = Vector{UInt8}(timezone)
        dpb = vcat(dpb, [isc_dpb_session_time_zone, UInt8(length(timezone))], timezone)
    end

    if length(wp.auth_data) != 0
        specific_auth_data = Vector{UInt8}(bytes2hex(wp.auth_data))
        dpb = vcat(dpb, [UInt8(isc_dpb_specific_auth_data), UInt8(length(specific_auth_data))], specific_auth_data)
    end

    if wp.timezone != ""
        tzname_bytes = Vector{UInt8}(wp.timezone)
        dpb = vcat(dpb, [isc_dpb_session_time_zone, length(tzname_bytes)], tzname_bytes)
    end

    pack_uint32(wp, op_attach)
    pack_uint32(wp, 0)
    pack_string(wp, db_name)
    pack_bytes(wp, dpb)
    send_packets(wp)
end

function _op_cont_auth(wp::WireProtocol, auth_data::Vector{UInt8}, auth_plugin_name::String, keys::String)
    DEBUG_OUTPUT("_op_cont_auth")
    pack_uint32(wp, op_cont_auth)
    pack_string(wp, bytes2hex(auth_data))
    pack_string(wp, auth_plugin_name)
    pack_string(wp, "Srp256,Srp")
    pack_string(wp, keys)
    send_packets(wp)
end

function _op_crypt(wp::WireProtocol, plugin)
    DEBUG_OUTPUT("_op_crypt")
    pack_uint32(wp, op_crypt)
    pack_string(wp, plugin)
    pack_string(wp, "Symmetric")
    send_packets(wp)
end

function _op_drop_database(wp::WireProtocol)
    DEBUG_OUTPUT("_op_drop_database")
    pack_uint32(wp, op_drop_database)
    pack_uint32(wp, wp.db_handle)
    send_packets(wp)
end

function _op_transaction(wp::WireProtocol, tpb::Vector{UInt8})
    DEBUG_OUTPUT("_op_transaction")
    pack_uint32(wp, op_transaction)
    pack_uint32(wp, wp.db_handle)
    pack_bytes(wp, tpb)
    send_packets(wp)
end

function _op_commit(wp::WireProtocol, trans_handle::Int32)
    DEBUG_OUTPUT("_op_commit")
    pack_uint32(wp, op_commit)
    pack_uint32(wp, trans_handle)
    send_packets(wp)
end

function _op_commit_retaining(wp::WireProtocol, trans_handle::Int32)
    DEBUG_OUTPUT("_op_commit_retaining")
    pack_uint32(wp, op_commit_retaining)
    pack_uint32(wp, trans_handle)
    send_packets(wp)
end

function _op_rollback(wp::WireProtocol, trans_handle::Int32)
    DEBUG_OUTPUT("_op_rollback")
    pack_uint32(wp, op_rollback)
    pack_uint32(wp, trans_handle)
    send_packets(wp)
end

function _op_rollback_retaining(wp::WireProtocol)
    DEBUG_OUTPUT("_op_rollback_retaining")
    pack_uint32(wp, op_rollback_retaining)
    pack_uint32(wp, trans_handle)
    send_packets(wp)
end

function _op_allocate_statement(wp::WireProtocol)
    DEBUG_OUTPUT("_op_allocate_statement")
    pack_uint32(wp, op_allocate_statement)
    pack_uint32(wp, wp.db_handle)
    send_packets(wp)
end

function _op_info_transaction(wp::WireProtocol, trans_handle::Int32, b::Vector{UInt8})
    DEBUG_OUTPUT("_op_info_transaction")
    pack_uint32(wp, op_info_transaction)
    pack_uint32(wp, trans_handle)
    pack_uint32(wp, 0)
    pack_bytes(wp, b)
    pack_uint32(wp, BUFFER_LEN)
    send_packets(wp)
end

function _op_info_database(wp::WireProtocol, b::Vector{UInt8})
    DEBUG_OUTPUT("_op_info_database")
    pack_uint32(wp, op_info_database)
    pack_uint32(wp, wp.db_handle)
    pack_uint32(wp, 0)
    pack_bytes(wp, b)
    pack_uint32(wp, BUFFER_LEN)
    send_packets(wp)
end

function _op_free_statement(wp::WireProtocol, stmt_handle::Int32, mode::Int)
    DEBUG_OUTPUT("_op_free_statement")
    pack_uint32(wp, op_free_statement)
    pack_uint32(wp, stmt_handle)
    pack_uint32(wp, mode)
    send_packets(wp)
end

function _op_prepare_statement(wp::WireProtocol, trans_handle::Int32, stmt_handle::Int32, query::String)
    DEBUG_OUTPUT("_op_prepare_statement")
    bs::Vector{UInt8} = vcat([isc_info_sql_stmt_type], INFO_SQL_SELECT_DESCRIBE_VARS())
    pack_uint32(wp, op_prepare_statement)
    pack_uint32(wp, trans_handle)
    pack_uint32(wp, stmt_handle)
    pack_uint32(wp, 3)  # dialect = 3
    pack_string(wp, query)
    pack_bytes(wp, bs)
    pack_uint32(wp, BUFFER_LEN)
    send_packets(wp)
end

function _op_info_sql(wp::WireProtocol, stmt_handle::Int32, vars::Vector{UInt8})
    DEBUG_OUTPUT("_op_info_sql")
    pack_uint32(wp, op_info_sql)
    pack_uint32(wp, stmt_handle)
    pack_uint32(wp, 0)
    pack_bytes(wp, vars)
    pack_uint32(wp, BUFFER_LEN)
    send_packets(wp)
end

function _op_execute(wp::WireProtocol, stmt_handle::Int32, trans_handle::Int32, params)
    DEBUG_OUTPUT("_op_execute")
    pack_uint32(wp, op_execute)
    pack_uint32(wp, stmt_handle)
    pack_uint32(wp, trans_handle)

    if length(params) == 0
        pack_uint32(wp, 0)  # pack_bytes(wp, [])
        pack_uint32(wp, 0)
        pack_uint32(wp, 0)
    else
        blr, values = params_to_blr(wp, trans_handle, params)
        pack_bytes(wp, blr)
        pack_uint32(wp, 0)
        pack_uint32(wp, 1)
        append_bytes(wp, values)
    end

    if wp.protocol_version >= PROTOCOL_VERSION16
        # statetment timeout
        append_bytes(wp, int32_to_bytes(Int32(0)))
    end

    send_packets(wp)

end

function _op_execute2(wp::WireProtocol, stmt_handle::Int32, trans_handle::Int32, params, output_blr::Vector{UInt8})
    DEBUG_OUTPUT("_op_execute2")
    pack_uint32(wp, op_execute2)
    pack_uint32(wp, stmt_handle)
    pack_uint32(wp, trans_handle)

    if length(params) == 0
        pack_uint32(wp, 0)  # pack_bytes(wp, [])
        pack_uint32(wp, 0)
        pack_uint32(wp, 0)
    else
        blr, values = params_to_blr(wp, trans_handle, params)
        pack_bytes(wp, blr)
        pack_uint32(wp, 0)
        pack_uint32(wp, 1)
        append_bytes(wp, values)
    end

    pack_bytes(wp, output_blr)
    pack_uint32(wp, 0)

    if wp.protocol_version >= PROTOCOL_VERSION16
        # statetment timeout
        append_bytes(wp, int32_to_bytes(Int32(0)))
    end

    send_packets(wp)
end

function _op_exec_immediate(wp::WireProtocol, trans_handle::Int32, query::String)
    DEBUG_OUTPUT("_op_exec_immediate")
    pack_uint32(wp, op_exec_immediate)
    pack_uint32(wp, trans_handle)
    pack_uint32(wp, wp.db_handle)
    pack_uint32(wp, 3)      # dialect = 3
    pack_string(wp, query)
    pack_uint32(wp, 0)      # pack_bytes(wp, [])
    pack_uint32(wp, BUFFER_LEN)
    send_packets(wp)
end

function _op_fetch(wp::WireProtocol, stmt_handle::Int32, blr::Vector{UInt8})
    DEBUG_OUTPUT("_op_fetch")
    pack_uint32(wp, op_fetch)
    pack_uint32(wp, stmt_handle)
    pack_bytes(wp, blr)
    pack_uint32(wp, 0)
    pack_uint32(wp, 400)
    send_packets(wp)
end

function _op_fetch_response(wp::WireProtocol, stmt_handle::Int32, xsqlda::Vector{XSQLVAR})::Tuple{Vector{Vector{Any}}, Bool}
    DEBUG_OUTPUT("_op_fetch_response")
    op_code = bytes_to_bint32(recv_packets(wp, 4))
    while op_code == op_dummy
        op_code = bytes_to_bint32(recv_packets(wp, 4))
    end
    while op_code == op_response && wp.lazy_response_count > 0
        wp.lazy_response_count -= 1
        op_code = bytes_to_bint32(recv_packets(wp, 4))
    end

    if op_code != op_fetch_response
        if op_code == op_response
            parse_op_response(wp)
        end
        throw(DomainError("op_fetch_resonse:op_code=$(op_code)"))
    end

    status = bytes_to_bint32(recv_packets(wp, 4))
    count = bytes_to_bint32(recv_packets(wp, 4))
    rows::Vector{Vector{Any}} = []

    while count > 0
        n = div(length(xsqlda), 8)
        if length(xsqlda) % 8 != 0
            n += 1
        end
        null_indicator::BigInt = 0
        for b in reverse(recv_packets_alignment(wp, n))
            null_indicator <<= 8
            null_indicator += b
        end
        r = Vector{Any}(missing, length(xsqlda))
        for i in 1:length(xsqlda)
            x = xsqlda[i]
            if (null_indicator & (1 << (i-1))) == 0 # not null
                ln = io_length(x)
                if ln < 0
                    ln = bytes_to_bint32(recv_packets(wp, 4))
                end
                raw_value = recv_packets_alignment(wp, ln)
                r[i] = value(x, raw_value)
            end
        end
        rows = vcat(rows, [r])
        op_code = bytes_to_bint32(recv_packets(wp, 4))
        status = bytes_to_bint32(recv_packets(wp, 4))
        count = bytes_to_bint32(recv_packets(wp, 4))
    end

    rows, status != 100
end

function _op_detatch(wp::WireProtocol)
    DEBUG_OUTPUT("_op_detatch")
    pack_uint32(wp, op_detatch)
    pack_uint32(wp, wp.db_handle)
    send_packets(wp)
end

function _op_open_blob(wp::WireProtocol, blob_id::Vector{UInt8}, trans_handle::Int32)
    DEBUG_OUTPUT("_op_open_blob")
    pack_uint32(wp, op_open_blob)
    pack_uint32(wp, trans_handle)
    append_bytes(wp, blob_id)
    send_packets(wp)
end

function _op_create_blob2(wp::WireProtocol, trans_handle::Int32)
    DEBUG_OUTPUT("_op_create_blob2")
    pack_uint32(wp, op_create_blob2)
    pack_uint32(wp, 0)
    pack_uint32(wp, trans_handle)
    pack_uint32(wp, 0)
    pack_uint32(wp, 0)
    send_packets(wp)
end

function _op_get_segment(wp::WireProtocol, blob_handle::Int32)
    DEBUG_OUTPUT("_op_get_segment")
    pack_uint32(wp, op_get_segment)
    pack_uint32(wp, blob_handle)
    pack_uint32(wp, BUFFER_LEN)
    pack_uint32(wp, 0)
    send_packets(wp)
end

function _op_put_segment(wp::WireProtocol, blob_handle::Int32, seg_data::Vector{UInt8})
    DEBUG_OUTPUT("_op_put_segment")
    ln::UInt32 = length(seg_data)

    pack_uint32(wp, op_put_segment)
    pack_uint32(wp, blob_handle)
    pack_uint32(wp, ln)
    pack_uint32(wp, ln)
    append_bytes(wp, seg_data)
    padding = Vector{UInt8}([0, 0, 0])
    append_bytes(wp, padding[1:((4 - ln) & 3)])
    send_packets(wp)
end

function _op_batch_segments(wp::WireProtocol, blob_handle::Int32, seg_data::Vector{UInt8})
    DEBUG_OUTPUT("_op_batch_segment")
    ln = length(sqg_data)
    pack_uint(wp, op_batch_segments)
    pack_uint32(wp, blob_handle)
    pack_uint32(wp, ln+2)
    pack_uint32(wp, ln+2)
    pad_length = ((4-(ln+2)) & 3)
    padding = Vector{UInt8}([0, 0, 0])
    pack_bytes(wp, Vector{UInt8}([UInt8(ln & 255), UInt8(ln >> 8)]))  # little endian int 16
    pack_bytes(wp, seg_data)
    append_bytes(wp, padding[1:((4 - ln) & 3)])
    send_packets(wp)
end

function _op_close_blob(wp::WireProtocol, blob_handle::Int32)
    DEBUG_OUTPUT("_op_close_blob")
    pack_uint32(wp, op_close_blob)
    pack_uint32(wp, blob_handle)
    send_packets(wp)
end

function _op_response(wp::WireProtocol)::Tuple{Int32, Vector{UInt8}, Vector{UInt8}}
    DEBUG_OUTPUT("_op_response")
    op_code = bytes_to_bint32(recv_packets(wp, 4))
    while op_code == op_dummy
        op_code = bytes_to_bint32(recv_packets(wp, 4))
    end
    while op_code == op_response && wp.lazy_response_count > 0
        wp.lazy_response_count -= 1
        parse_op_response(wp)
        op_code = bytes_to_bint32(recv_packets(wp, 4))
    end
    if op_code == op_cont_auth
        throw(DomainError("Unauthorized"))
    elseif op_code != op_response
        throw(DomainError("op_resonse:op_code=$(op_code)"))
    end
    parse_op_response(wp)
end

function _op_sql_response(wp::WireProtocol, xsqlda::XSQLVAR)::Vector{Any}
    DEBUG_OUTPUT("_op_sql_response")
    op_code = bytes_to_bint32(recv_packets(wp, 4))
    while op_code == op_dummy
        op_code = bytes_to_bint32(recv_packets(wp, 4))
    end
    if op_code != op_sql_response
        if op_code == op_response
            parse_op_response(wp)
        end
        throw(DomainError("op_sql_response:op_code=$(op_code)"))
    end

    r::Vector{Any} = []
    count = bytes_to_bint32(recv_packets(wp, 4))
    if count == 0
        return r
    end

    n = div(length(xsqlvar), 8)
    if length(xsqlvar) % 8 != 0
        n += 1
    end
    null_indicator = 0
    for c in reverse(recv_packets_alignment(wp, n))
        null_indicator <<= 8
        null_indicator += 8
    end
    for i in 1:length(xsqlda)
        x = xsqlda[i]
        if (null_indicator & (1 << (i-1))) != 0
            append!(r, nothing)
        else
            ln = io_length(x)
            if ln < 0
                ln = bytes_to_bint32(recv_packets(wp, 4))
            end
            raw_value = recv_packets_alignment(wp, ln)
            append!(r, value(x, raw_value))
        end
    end
    r
end

function create_blob(wp::WireProtocol, b::Vector{UInt8}, trans_handle::Int32)::Vector{UInt8}
    buf = suspebd_buffer(wp)
    resume_buffer(wp, buf)

    _op_create_blob2(wp, trans_handle)
    blob_handle, blob_id, _ = op_response(wp)

    i = 1
    while i < length(b)
        _op_put_segent(wp, blob_handle, b[i:i+BLOB_SEGMENT_SIZE])
        h, oid, buf = op_response(wp)
        i += BLOB_SEGMENT_SIZE
    end
    op_close_blob(wp, blob_handle)
    h, oid, buf = op_response(wp)

    resule_buffer(buf)

    blob_id
end

function params_to_blr(wp::WireProtocol, trans_handle::Int32, params)::Tuple{Vector{UInt8}, Vector{UInt8}}
    ln = length(params) * 2
    blr_list::Vector{UInt8} = Vector{UInt8}([5, 2, 4, 0, UInt8(ln & 255), UInt8(ln >> 8)])
    values::Vector{UInt8} = []

    # NULL indicator
    null_indicator = 0
    for i in 1:length(params)
        if params[i] == nothing
            null_indicator |= (1 << i)
        end
    end
    n = div(length(params), 8)
    if length(params) % 8 != 0
        n += 1
    end
    if n % 4 != 0   # padding
        n += (4 - n) % 4
    end
    for i in 1:n
        append!(values, null_indicator & 255)
        null_indicator >>= 8
    end

    for i in 1:length(params)
        p = params[i]
        if typeof(p) == String
            p = Vector{UInt8}(params[i])
        end
        v, blr = if typeof(p) == Vector{UInt8}
            if length(p) < MAX_CHAR_LENGTH
                bytes_to_blr(p)
            else
                create_blob(wp, p, trans_handle), UInt8[9, 0]
            end
        else
            to_blr(p)
        end
        append!(values, v)
        append!(blr_list, blr)
        append!(blr_list, 7)
        append!(blr_list, 0)
    end

    blr_list = vcat(blr_list, Vector{UInt8}([255, 76]))     # [blr_end, blr_eoc]
    blr_list, values
end


text:
using Cuddle
using Base.Test

println("Massive cuddle sesh initiated. Sweet.")

ext: jl


text:
include("JuliaNature.jl")
f(x1, x2) = (x1 + 2 * x2 - 7)^2 + (2 * x1 + x2 - 5)^2
model = JuliaNature.initialize(25, 2, f, [0, 1.0], [-10.0, 10.0], [0.0, 2.0], [0, 1.0], 0)
println(JuliaNature.optimize(model, 10000))

max_stars_repo_name: JuliaNature/JuliaNature

max_issues_repo_name: JuliaNature/JuliaNature

max_forks_repo_name: JuliaNature/JuliaNature


text:
export BatchedSampler

mutable struct SimpleIterativeSampler{BN, S, Sc, Sv, H, Pd} <: AbstractIterativeSampler
    bnet::BN

    sampler::S
    sampler_cache::Sc
    samples::Sv

    hilb::H

    parallel_cache::Pd
end

function SimpleIterativeSampler(net,
                                sampl,
                                hilb;
                                batch_sz=2^4,
                                par_type=automatic_parallel_type())
    if net isa CachedNet
        throw("Only takes standard network.")
    end

    if hilb isa AbstractOperator
        hilb = basis(hilb)
    end

    par_cache      = parallel_execution_cache(par_type)

    bnet           = cached(net, batch_sz)
    v              = state(hilb, bnet)
    sampler_cache  = cache(sampl, v, hilb, bnet, par_cache)

    ch_len         = chain_length(sampl, sampler_cache)

    samples        = state(hilb, bnet, ch_len)

    sis = SimpleIterativeSampler(bnet,
            sampl, sampler_cache, samples,
            hilb,
            par_cache)

    return sis
end

function sample!(is::SimpleIterativeSampler)
    # Sample the new configurations
    _sample_state!(is)

    return is.samples
end

samples(is::SimpleIterativeSampler) = is.samples
chain_length(is::SimpleIterativeSampler) = chain_length(is.sampler, is.sampler_cache)

Base.show(io::IO, is::SimpleIterativeSampler) = print(io,
    "SimpleIterativeSampler for :"*
    "\n\tnet\t\t: $(is.bnet)"*
    "\n\tsampler\t: $(is.sampler)")


text:
include("../common/PWGrid_v03.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")
include("gen_dr.jl")
include("init_pot_harm_3d.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_sd.jl")
include("KS_solve_Emin_cg.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

function test_main( ecutwfc_Ry::Float64 )

    const LatVecs = 6.0*diagm( ones(3) )

    pw = PWGrid( 0.5*ecutwfc_Ry, LatVecs )

    const Ω  = pw.Ω
    const r  = pw.r
    const G  = pw.gvec.G
    const G2 = pw.gvec.G2
    const Ns = pw.Ns
    const Npoints = prod(Ns)
    const Ngwx = pw.gvecw.Ngwx

    @printf("ecutwfc (Ha) = %f\n", ecutwfc_Ry*0.5)
    @printf("Ns = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)

    const actual = Npoints/Ngwx
    const theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    #
    # Generate array of distances
    #
    center = sum(LatVecs,2)/2
    dr = gen_dr( r, center )
    #
    # Setup potential
    #
    V_ionic = init_pot_harm_3d( pw, dr )
    print("sum(Vpot)*Ω/Npoints = $(sum(V_ionic)*Ω/Npoints)\n");
    #
    const Nstates = 4
    Focc = 2.0*ones(Nstates)
    #
    #psi, Energies, Potentials = KS_solve_Emin_sd( pw, V_ionic, Focc, Nstates, NiterMax=10 )
    #psi, Energies, Potentials = KS_solve_Emin_cg( pw, V_ionic, Focc, Nstates,
    #                            NiterMax=1000, Potentials0=Potentials, psi0=psi )

    psi, Energies, Potentials = KS_solve_Emin_cg( pw, V_ionic, Focc, Nstates, NiterMax=1000 )

    #
    print_Energies(Energies)

    Y = ortho_gram_schmidt(psi)
    mu = Y' * op_H( pw, Potentials, Y )
    evals, evecs = eig(mu)
    Psi = Y*evecs
    for st = 1:Nstates
        @printf("=== State # %d, Energy = %f ===\n", st, real(evals[st]))
    end
end

@time test_main( 40.0 )


text:
# the purpose of this file is to find a winner for the fastest histogram
# It was found that (unsurprisingly) the multithreaded version is faster
# actually the story is a bit murky now and it's best to wait for v0.7.0
# I seem to find type instability issues with the threadedcount code
# but it disappeared in v0.7.0, so waiting might be a good strategy
tic()
using Base.Threads
using Plots
using BenchmarkTools, DataFrames
import SortingAlgorithms: uint_mapping, RADIX_MASK, RADIX_SIZE, load_bits

# multi-threaded histogram count
# there seems to be some instability
function threadedcount(bits::Vector{T}) where T <: Unsigned
    iter = sizeof(T)÷2::Int
    hist = zeros(UInt32, 65536, iter, nthreads())
    l = length(bits)::Int
    @threads for j = 1:l
        # for i = 0:iter-1
        #     # idx = (bits[j] >> (i << 4)) & T(0xffff)
        #     # idx = bits[j]
        #     # @inbounds hist[idx, i+1, threadid()] += 1
        # end
    end
    nt = Threads.nthreads()
    @threads for j in 1:iter
        for i = 2:nt
            # @inbounds hist[:, j, 1] .+= hist[:, j, i]
        end
    end
    hist[:,:,1]
end

@code_warntype threadedcount()

# single threaded histogram count from SortingAlgorithms.jl
function hist_sortingalgorithms(bits::Vector{T}) where T<:Unsigned
    iters = sizeof(T)÷2
    bin = zeros(UInt32, 65536, iters)

    # Histogram for each element, radix
    for v = bits
        for j = 1:iters
            idx = Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK) + 1
            @inbounds bin[idx,j] += 1
        end
    end
end

function count64(bits::Vector{T}) where T<:Unsigned
    hist = zeros(UInt32, 1<<16, 4)
    for x1 in bits
        for i in 0:sizeof(T)÷2-1
            @inbounds hist[1+Int((x1 >> (i << 4)) & 0xffff), i+1] += 1
        end
    end
    hist
end

function count16(bits::Vector{T}) where T<:Unsigned
    hist = zeros(UInt32, 1<<16, 4)
    hist_hot = zeros(UInt16, 1<<16)
    @inbounds for i in 0:sizeof(T)÷2-1
        for x1 in bits
            idx = 1+ Int((x1 >> (i << 4)) & 0xffff)
            hist_hot[idx] += UInt16(1)
            if hist_hot[idx] == 0
                hist[idx, i+1] += (1<<16)
            end
        end
        for idx in 1:(1<<16)
            hist[idx,i+1] += hist_hot[idx]
        end
        fill!(hist_hot, zero(UInt16))
    end
    hist
end

function count8(bits::Vector{T}) where T<:Unsigned
    hist = zeros(UInt32, 1<<16, 4)
    hist_hot = zeros(UInt8, 1<<16)
    @inbounds for i in 0:sizeof(T)÷2-1
        for x1 in bits
            idx = 1+ Int((x1 >> (i << 4)) & 0xffff)
            hist_hot[idx] += UInt8(1)
            if hist_hot[idx] == 0
                hist[idx, i+1] += (1<<8)
            end
        end
        for idx in 1:(1<<16)
            hist[idx,i+1] += hist_hot[idx]
        end
        fill!(hist_hot, zero(UInt8))
    end
    hist
end

function count8t(bits::Vector{T}) where T<:Unsigned
    hist = zeros(UInt32, 1<<16, sizeof(T)÷2)
    hist_hot = zeros(UInt8, 1<<16, sizeof(T)÷2)
    @threads for i in 0:sizeof(T)÷2-1
        @inbounds for x1 in bits
            idx = 1+ Int((x1 >> (i << 4)) & 0xffff)
            hist_hot[idx,i+1] += UInt8(1)
            if hist_hot[idx,i+1] == 0
                hist[idx, i+1] += (1<<8)
            end
        end
        @inbounds for idx in 1:(1<<16)
            hist[idx,i+1] += hist_hot[idx, i+1]
        end
    end
    hist
end

function count8x(bits::Vector{T}) where T<:Unsigned
    hists = [zeros(UInt32, 1<<16) for i in 1:sizeof(T)÷2]
    hists_hot = [zeros(UInt8, 1<<16) for i in 1:sizeof(T)÷2]
    @threads for i in 0:sizeof(T)÷2-1
        hist = hists[i+1]
        hist_hot = hists_hot[i+1]
        i4 = UInt8(i<<4)
        @inbounds for x1 in bits
            idx = 1+ Int((x1 >> i4) & 0xffff)
            hist_hot[idx] += UInt8(1)
            if hist_hot[idx] == 0
                hist[idx] += (1<<8)
            end
        end
        @inbounds for idx in 1:(1<<16)
            hist[idx] += hist_hot[idx]
        end
    end
    hcat(hists...)
end

using BenchmarkTools
N = 100_000_000
K = 100

function test_histogram(N,K)
    @time samplespace = "id".*dec.(1:N÷K,10);
    @time svec = rand(load_bits.(UInt32, samplespace), N);

    naiveres = @belapsed threadedcount($svec)
    sortingalgorithms = @belapsed hist_sortingalgorithms($svec)
    count64res = @belapsed count64($svec)
    count16res = @belapsed count16($svec)
    count8res = @belapsed count8($svec)
    count8tres = @belapsed count8t($svec)
    count8xres = @belapsed count8x($svec)

    first4res = DataFrame(
        alg = ["orig",    "SortingAlgorithms.jl", "count64",  "count16",  "count8",
            "count8t",  "count8x"],
        timing = [naiveres,  sortingalgorithms,      count64res, count16res, count8res,
            count8tres,  count8xres])

    b1 = bar(first4res[:alg], first4res[:timing],
        label="time (sec)", title = "sort first 4 bytes - $(N/1_000_000)m")

    @time svec = rand(load_bits.(UInt, samplespace, 4), N);
    gc()
    naiveres = @belapsed threadedcount($svec)
    sortingalgorithms = @belapsed hist_sortingalgorithms($svec)
    count64res = @belapsed count64($svec)
    count16res = @belapsed count16($svec)
    count8res = @belapsed count8($svec)
    count8tres = @belapsed count8t($svec)
    count8xres = @belapsed count8x($svec)

    last8res = DataFrame(
        alg = ["orig",    "SortingAlgorithms.jl", "count64",  "count16",  "count8",
            "count8t",  "count8x"],
        timing = [naiveres,  sortingalgorithms,      count64res, count16res, count8res,
            count8tres,  count8xres])

    b2 = bar(last8res[:alg], last8res[:timing],
        label="time (sec)", title = "sort last 8 bytes - $(N/1_000_000)m")

    plot(b1,b2)
    savefig("compare_rbench_count_uint.png")

    first4res, last8res, b1, b2
end

test_histogram(1_000_000, 100)

res = test_histogram(100_000_000, 100)
toc()

N = 1000
K = 100
import SortingLab: load_bits
@time samplespace = "id".*dec.(1:N÷K,10);
@time svec = rand(load_bits.(UInt32, samplespace), N);
@code_warntype threadedcount(svec)

# @btime threadedcount($svec);
# @btime threadedcount($svec, UInt32); # a bit slower but uses only half the memory
# svec = rand(rand(UInt, 10_000_000), 1_000_000_000);
# d = @elapsed threadedcount(svec)
# a = @elapsed count64(svec)
# b = @elapsed count16(svec)
# c1 = @elapsed count8(svec)
# c = @elapsed count8t(svec)
# f = @elapsed count8x(svec)

# savefig("Single vs Multithreaded countall.png")

# using Plots
# using BenchmarkTools
# svec = rand(rand(UInt(1):UInt(10_000_000), 10_000_000), 1_000_000_000);
# d = @elapsed threadedcount(svec)
# a = @elapsed count64(svec)
# b = @elapsed count16(svec)
# c1 = @elapsed count8(svec)
# c = @elapsed count8t(svec)
# f = @elapsed count8x(svec)

# bar(["orig","count64", "count16","count8", "count8t", "count8x"], [d, a, b,c1, c, f], label="time (sec)")
# savefig("Single vs Multithreaded count.png")

# bar(["orig","count16","count8", "count8t", "count8x"], [d, b,c1, c, f], label="time (sec)")
# savefig("Single vs Multithreaded count.png")

# function abc(bits::Vector{T}) where T <: Unsigned
#     iter = sizeof(T)÷2::Int
#     # hist = zeros(UInt32, 65536, iter, nthreads())
#     # l = length(bits)::Int
#     # @threads for j = 1:l
#     #     # for i = 0:iter-1
#     #     #     # idx = (bits[j] >> (i << 4)) & T(0xffff)
#     #     #     # idx = bits[j]
#     #     #     # @inbounds hist[idx, i+1, threadid()] += 1
#     #     # end
#     # end
#     nt = Threads.nthreads()
#     @threads for k in 1:iter
#         for i = 2:nt
#             # @inbounds hist[:, j, 1] .+= hist[:, j, i]
#         end
#     end
#     # hist[:,:,1]
#     1
# end

# @code_warntype abc(svec)

# function def()
#     @threads for j = 1:8
#     end

#     @threads for k = 1:8
#     end

#     nothing
# end

# @code_warntype def()


text:
"""
  Abstract type for all shock capturing methods.  This type determines what
  kind of dissipation to add when there is a shock present.  Ideally, any
  shock capturing scheme should be able to be paired with any shock sensor.

  New shock capturing scheme should generally not subtype this type directly,
  instead they should subtype either [`AbstractVolumeShockCapturing`](@ref)
  or [`AbstractFaceShockCapturing`](@ref).

  Optional methods:

  ```
    getShockSensor
    getShockSensor
    getDiffusion
  ```

  `getShockSensor` and `setShockSensor only need to be extended if the
  `AbstractShockCapturing` does not store the sensor in the `.sensor` field.

  `getDiffusion` only needs to be extended if the `AbstractShockCapturing` does
  not store the sensor in the `.diffusion` field.

"""
abstract type AbstractShockCapturing end

"""
  Returns the shock sensor.  This function is not required to be type-stable

  **Inputs**

   * obj: an `AbstractShockCapturing`

  **Outputs**

   * sensor: an `AbstractShockSensor`
"""
function getShockSensor(obj::AbstractShockCapturing)

  return obj.sensor
end

"""
  Returns the `AbstractDiffusion` object.  This function is not required to
  be type stable.

  **Inputs**

   * obj: an `AbstractShockCapturing`

  **Outputs**

   * diffusion: an `AbstractDiffusion`
"""
function getDiffusion(obj::AbstractShockCapturing)

  return obj.diffusion
end

"""
  Counterpart of `getShockSensor`, sets a new shock sensor to be used with
  the existing `AbstractShockCapturing`.

  **Inputs**

   * obj: `AbstractShockCapturing`
   * sensor: the new `AbstractShockSensor`
"""
function setShockSensor(obj::AbstractShockCapturing, sensor::AbstractShockSensor)
  obj.sensor = sensor
  obj.diffusion = ShockDiffusion(obj.diffusion, sensor)
end

#------------------------------------------------------------------------------
# Volume shock capturing

"""
  Abstract type for shock capturing methods that only have volume terms.
  These scheme are substantially simpler and do not require constructing
  a reduced mesh.  The interface for this type has not been finalized

  This type must implement functions:

  ```
    calcShockCapturing
    calcShockCapturing_diff
  ```

  and may optionally implement

  ```
    getSparsityPattern
  ```
"""
abstract type AbstractVolumeShockCapturing <: AbstractShockCapturing end

"""
  This function should update the residual with the shock capturing term

  **Inputs**

   * mesh: mesh object
   * sbp: SBP operator
   * eqn: EulerData
   * opts: options dictonary
   * sensor: [`AbstractShockSensor`](@ref)
   * capture: the shock capturing scheme
"""
function calcShockCapturing(mesh::AbstractMesh, sbp::AbstractOperator,
                            eqn::EulerData, opts,
                            sensor::AbstractShockSensor,
                            capture::AbstractVolumeShockCapturing)

  error("abstract fallback for calcShockCapturing() called. Did you forget to extend it with a new method for your AbstractFaceShockCapturing?")

end

"""
  This function should compute the Jacobian of the shock capturing term

  **Inputs**

   * mesh: mesh object
   * sbp: SBP operator
   * eqn: EulerData
   * opts: options dictionary
   * sensor: an [`AbstractShockSensor`](@ref)
   * capture: the shock capturing object
   * assem: an [`AssembleElementData`](@ref)

  Note that many `capture` types require the same sensor to be passed into
  their constructor as to `calcShockCapturing`.

"""
function calcShockCapturing_diff(mesh::AbstractMesh, sbp::AbstractOperator,
                                 eqn::EulerData, opts,
                                 sensor::AbstractShockSensor,
                                 capture::AbstractVolumeShockCapturing,
                                 assem::AssembleElementData)

  error("abstract fallback for calcShockCapturing_diff() called. Did you forget to extend it with a new method for your AbstractFaceShockCapturing?")
end

"""
  This function returns an enum specifying the sparsity pattern of the
  Jacobian required by this shock capturing scheme.  By default,
  `AbstractVolumeShockCapturing` has an `INVISCID` sparsity pattern and
  `AbstractFaceShockCApturing` has a `VISCOUSTIGHT` sparsity pattern.

  Shock capturing scheme can introduce a specailized method if the defaults
  do not work for them.

  **Inputs**

   * capture: an `AbstractShockCapturing`

  **Outputs**

   * enum specifying sparsity pattern
"""
function getSparsityPattern(capture::AbstractVolumeShockCapturing)

  return INVISCID
end

#------------------------------------------------------------------------------
# AbstractFaceShockCapturing

"""
  Abstract type for shock capturing methods that do face integrals (these
  scheme may do volume integrals as well).  These schemes require constructing
  a reduced mesh for the elements that have shocks in them.

  These types must implement these functions:
  ```
    allocateArrays
    calcShockCapturing
    calcShockCapturing_diff
  ```

"""
abstract type AbstractFaceShockCapturing <: AbstractShockCapturing end

"""
  This function is called every time the `shockmesh` is updated, and performs
  any setup work that depends on the `shockmesh`

  **Inputs**

   * capture: [`LDGShockCapturing`](@ref)
   * mesh
   * shockmesh: a `ShockedElements` object, fully initialized

"""
function allocateArrays(capture::AbstractFaceShockCapturing,
                        mesh::AbstractMesh, shockmesh::ShockedElements)

  error("abstract fallback for allocateArrays() called. Did you forget to extend it with a new method for your AbstractFaceShockCapturing?")
end

"""
  Update `eqn.res` with the shock capturing term

  **Inputs**

   * mesh
   * sbp
   * eqn: `eqn.res` should be updated
   * opts
   * capture: the shock capturing scheme to be used (this argument must
              be specialized)
   * shockmesh

"""
function calcShockCapturing(mesh::AbstractMesh, sbp::AbstractOperator,
                           eqn::EulerData, opts,
                           capture::AbstractFaceShockCapturing,
                           shockmesh::ShockedElements)

  error("abstract fallback for calcShockCapturing() called. Did you forget to extend it with a new method for your AbstractFaceShockCapturing?")
end

"""
  Assemble the Jacobian of the shock capturing term into the Jacobian

  **Inputs**

   * mesh
   * sbp
   * eqn:
   * opts
   * capture: the shock capturing scheme to be used (this argument must
              be specialized)
   * shockmesh: a `ShockedElements`
   * assem:  an `AssembleElementData`

"""
function calcShockCapturing_diff(mesh::AbstractMesh, sbp::AbstractOperator,
                             eqn::EulerData, opts,
                             capture::AbstractFaceShockCapturing,
                             shockmesh::ShockedElements,
                             assem::AssembleElementData)

  error("abstract fallback for calcShockCapturing_diff() called. Did you forget to extend it with a new method for your AbstractFaceShockCapturing?")
end

function getSparsityPattern(capture::AbstractFaceShockCapturing)

  return VISCOUSTIGHT
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

import Base.Threads.@threads
using TinyBnB, CPLEXW, Printf

instance_name = ARGS[1]
output_filename = ARGS[2]
node_limit = parse(Int, ARGS[3])

mip = open_mip(instance_name)
n_vars = CPXgetnumcols(mip.cplex_env[1], mip.cplex_lp[1])

pseudocost_count_up = [0 for i in 1:n_vars]
pseudocost_count_down = [0 for i in 1:n_vars]
pseudocost_sum_up = [0. for i in 1:n_vars]
pseudocost_sum_down = [0. for i in 1:n_vars]

function full_strong_branching_track(node::Node, progress::Progress)::TinyBnB.Variable
    N = length(node.fractional_variables)
    scores = Array{Float64}(undef, N)
    rates_up = Array{Float64}(undef, N)
    rates_down = Array{Float64}(undef, N)

    @threads for v in 1:N
        fix_vars!(node.mip, node.branch_variables, node.branch_values)
        obj_up, obj_down = TinyBnB.probe(node.mip, node.fractional_variables[v])
        unfix_vars!(node.mip, node.branch_variables)
        delta_up  = obj_up - node.obj
        delta_down = obj_down - node.obj
        frac_up = ceil(node.fractional_values[v]) - node.fractional_values[v]
        frac_down = node.fractional_values[v] - floor(node.fractional_values[v])
        rates_up[v] = delta_up / frac_up
        rates_down[v] = delta_down / frac_down
        scores[v] = delta_up * delta_down
    end

    max_score, max_offset = findmax(scores)
    selected_var = node.fractional_variables[max_offset]

    if abs(rates_up[max_offset]) < 1e6
        pseudocost_count_up[selected_var.index] += 1
        pseudocost_sum_up[selected_var.index] += rates_up[max_offset]
    end

    if abs(rates_down[max_offset]) < 1e6
        pseudocost_count_down[selected_var.index] += 1
        pseudocost_sum_down[selected_var.index] += rates_down[max_offset]
    end

    return selected_var
end

branch_and_bound(mip,
                 node_limit = node_limit,
                 branch_rule = full_strong_branching_track,
                 node_rule = best_bound,
                 print_interval = 100)

priority = [(pseudocost_count_up[v] == 0 || pseudocost_count_down[v] == 0) ? 0 :
              (pseudocost_sum_up[v] / pseudocost_count_up[v]) *
              (pseudocost_sum_down[v] / pseudocost_count_down[v])
             for v in 1:n_vars];

open(output_filename, "w") do file
    for var in mip.binary_variables
        write(file, @sprintf("%s,%.0f\n", name(mip, var), priority[var.index]))
    end
end

max_stars_repo_name: bknueven/MIPLearn

max_issues_repo_name: bknueven/MIPLearn

max_forks_repo_name: bknueven/MIPLearn


text:
"""

Help composite type as used to construct the cmdline.

$(SIGNATURES)

# Extended help

Primarily used for testing of StanBase. `Help()` defaults to `Help(:sample)`.

Example: `Help(:optimize)`
"""
struct Help
  help::Symbol
end
Help() = Help(:sample)


text:
using Primes

function count_prime_sums(limit::Int)
    p_list = primes(limit)
    counts = Dict(0=>1)

    for p in p_list
        for i in p:limit
            counts[i] = get(counts, i, 0) + get(counts, i - p, 0)
        end
    end

    for i in 2:limit
        if counts[i] > 5000
            return i
        end
    end
end

println(count_prime_sums(100))

max_stars_repo_name: HarrisonGreen/Project-Euler-Solutions

max_issues_repo_name: HarrisonGreen/Project-Euler-Solutions

max_forks_repo_name: HarrisonGreen/Project-Euler-Solutions


text:
mutable struct CuNFFTPlan{T,D} <: AbstractNFFTPlan{T,D,1}
  N::NTuple{D,Int64}
  NOut::NTuple{1,Int64}
  M::Int64
  x::Matrix{T}
  n::NTuple{D,Int64}
  dims::UnitRange{Int64}
  params::NFFTParams{T}
  forwardFFT::CUDA.CUFFT.cCuFFTPlan{Complex{T},-1,true,D}
  backwardFFT::CUDA.CUFFT.cCuFFTPlan{Complex{T},1,true,D}
  tmpVec::CuArray{Complex{T},D}
  tmpVecHat::CuArray{Complex{T},D}
  deconvolveIdx::CuArray{Int64,1}
  windowLinInterp::Vector{T}
  windowHatInvLUT::CuArray{Complex{T}} # ::Vector{Vector{T}}
  B::CuSparseMatrixCSC{Complex{T}} # ::SparseMatrixCSC{T,Int64}
end

function AbstractNFFTs.plan_nfft(::Type{<:CuArray}, x::Matrix{T}, N::NTuple{D,Int}, rest...;
  timing::Union{Nothing,TimingStats} = nothing, kargs...) where {T,D}
  t = @elapsed begin
    p = CuNFFTPlan(x, N, rest...; kargs...)
  end
  if timing != nothing
    timing.pre = t
  end
  return p
end

function CuNFFTPlan(x::Matrix{T}, N::NTuple{D,Int}; dims::Union{Integer,UnitRange{Int64}}=1:D,
                 fftflags=nothing, kwargs...) where {T,D}

    if dims != 1:D
      error("CuNFFT does not work along directions right now!")
    end

    params, N, NOut, M, n, dims_ = NFFT.initParams(x, N, dims; kwargs...)
    params.storeDeconvolutionIdx = true # CuNFFT only works this way
    params.precompute = NFFT.FULL # CuNFFT only works this way

    tmpVec = CuArray{Complex{T},D}(undef, n)

    #fftflags_ = (fftflags != nothing) ? (flags=fftflags,) : NamedTuple()
    #FP = plan_fft!(tmpVec, dims_; fftflags_...)
    #BP = plan_bfft!(tmpVec, dims_; fftflags_...)
    FP = plan_fft!(tmpVec, dims_)
    BP = plan_bfft!(tmpVec, dims_)

    windowLinInterp, windowHatInvLUT, deconvolveIdx, B = NFFT.precomputation(x, N[dims_], n[dims_], params)

    U = params.storeDeconvolutionIdx ? N : ntuple(d->0,D)
    tmpVecHat = CuArray{Complex{T},D}(undef, U)

    deconvIdx = CuArray(deconvolveIdx)
    winHatInvLUT = CuArray(Complex{T}.(windowHatInvLUT[1]))
    B_ = CuSparseMatrixCSC(Complex{T}.(B))

    CuNFFTPlan{T,D}(N, NOut, M, x, n, dims_, params, FP, BP, tmpVec, tmpVecHat,
               deconvIdx, windowLinInterp, winHatInvLUT, B_)
end

AbstractNFFTs.size_in(p::CuNFFTPlan) = p.N
AbstractNFFTs.size_out(p::CuNFFTPlan) = p.NOut

"""  in-place NFFT on the GPU"""
function LinearAlgebra.mul!(fHat::CuArray, p::CuNFFTPlan{T,D}, f::CuArray;
                          verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {T,D}
  NFFT.consistencyCheck(p, f, fHat)

  # deconvolve
  fill!(p.tmpVec, zero(Complex{T}))
  t1 = @elapsed begin
    p.tmpVecHat[:] .= vec(f).*p.windowHatInvLUT
    p.tmpVec[p.deconvolveIdx] = p.tmpVecHat
  end

  # FFT
  t2 = @elapsed p.forwardFFT * p.tmpVec

  # convolution
  t3 = @elapsed mul!(fHat, transpose(p.B), vec(p.tmpVec))

  if verbose
    @info "Timing: deconv=$t1 fft=$t2 conv=$t3"
  end
  if timing != nothing
    timing.conv = t3
    timing.fft = t2
    timing.deconv = t1
  end

  return fHat
end

"""  in-place adjoint NFFT on the GPU"""
function LinearAlgebra.mul!(f::CuArray, pl::Adjoint{Complex{T},<:CuNFFTPlan{T,D}}, fHat::CuArray;
                       verbose=false, timing::Union{Nothing,TimingStats} = nothing) where {T,D}
  p = pl.parent

  NFFT.consistencyCheck(p, f, fHat)

  # adjoint convolution
  t1 = @elapsed mul!(vec(p.tmpVec), p.B, fHat)

  # FFT
  t2 = @elapsed p.backwardFFT * p.tmpVec

  # adjoint deconvolve
  t3 = @elapsed begin
    p.tmpVecHat[:] = p.tmpVec[p.deconvolveIdx]
    f[:] .= vec(p.tmpVecHat) .* p.windowHatInvLUT
  end

  if verbose
    @info "Timing: conv=$t1 fft=$t2 deconv=$t3"
  end
  if timing != nothing
    timing.conv_adjoint = t1
    timing.fft_adjoint = t2
    timing.deconv_adjoint = t3
  end

  return f
end


text:
using Dash
using DashCoreComponents
using DashHtmlComponents
using PlotlyJS
include("IntegralUtils.jl")
include("GraphingUtils.jl")
using .IntegralUtils
using .GraphingUtils

mathjax = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js?config=TeX-MML-AM_CHTML"
app = dash(external_stylesheets=["https://codepen.io/chriddyp/pen/bWLwgP.css"],
            assets_folder="assets",
            external_scripts=[mathjax])

options_ = ["S: r(u, v)", "S: z = f(x, y), y = g(x)"]

app.layout = html_div() do
    html_div(children = [
    html_h1("\$\$\\text{Surface integrals calculator}\$\$", style=Dict("textAlign" => "center", "margin" => "3%")),
    html_div(id="header",
        children=[
        html_h3(id="header_h3", "\$\$ \\text{Hello, let the fun begin!}\$\$")
                ]
            ),
    html_hr(style=Dict("borderTop" => "1px dashed #d69c2f")),
    html_div(id="surface", children=[
        html_label(id="surface_label", "\$\$\\text{Choose the way of defining the surface:}\$\$"),
        dcc_dropdown(
            id = "dropdown",
            options = [(label = i, value = i) for i in options_],
            value = options_[1],
            ),
    ]),
    html_hr(style=Dict("borderTop" => "1px dashed #d69c2f")),
    html_div(id="main", children=[
        html_div(id="inputs", children=[
            html_label(id="vector_field_l", "\$\$\\text{Enter your vector field formula}\$\$"),
            html_div(id="field", children=[
            html_label(id="mFx", "\$\$\\vec{F}\\cdot\\hat{i}\$\$"),
            dcc_input(id="Fx", type="text", value="x"),
            html_label(id="mFy", "\$\$\\vec{F}\\cdot\\hat{j}\$\$"),
            dcc_input(id="Fy", type="text", value="y"),
            html_label(id="mFz", "\$\$\\vec{F}\\cdot\\hat{k}\$\$"),
            dcc_input(id="Fz", type="text", value="z"),
            ]
        ),
        html_label(id="method_dropdown", "\$\$\\text{Choose method of numeric integration:}\$\$"),
        dcc_dropdown(
            id = "method",
            options = [(label = i, value = i) for i in ("Simpson", "Monte Carlo")],
            value = "Simpson"
            ),
        html_label(id="accurary", "\$\$\\text{Choose the accuracy of numeric integration:}\$\$"),
        dcc_slider(
            id = "integral_accuracy",
            min = 10,
            max = 100,
            marks = Dict([Symbol(v) => Symbol(v) for v in 10:10:100]),
            value = 10,
            step = 10,
            ),
            html_div(id="parametric_bounds", children=[
            html_label(id="r1lab", "\$\$\\vec{r}(u, v)\\cdot\\hat{i} =\$\$"),
            dcc_input(id="r1", type="text", value="sqrt(1/4 + u^2) * cos(v)"),
            html_label(id="r2lab", "\$\$\\vec{r}(u, v)\\cdot\\hat{j} =\$\$"),
            dcc_input(id="r2", type="text", value="sqrt(1/4 + u^2) * sin(v)"),
            html_label(id="r3lab", "\$\$\\vec{r}(u, v)\\cdot\\hat{k} =\$\$"),
            dcc_input(id="r3", type="text", value="u"),
            html_label(id="u_range1", "\$\$u_{min} =\$\$"),
            dcc_input(id="u_range1i", type="text", value="-1"),
            html_label(id="u_range2", "\$\$u_{max} =\$\$"),
            dcc_input(id="u_range2i", type="text", value="1"),
            html_label(id="v_range1", "\$\$\\phi(u) =\$\$"),
            dcc_input(id="v_range1i", type="text", value="0"),
            html_label(id="v_range2", "\$\$\\theta(u) =\$\$"),
            dcc_input(id="v_range2i", type="text", value="2pi")
            ], style=Dict("display" => "grid")
            ),
        html_div(id="fxy_bounds", children=[
            html_label(id="x_range1", "\$\$x_{\\text{min}} =\$\$"),
            dcc_input(id="x_range1i", type="text", value="-2"),
            html_label(id="x_range2", "\$\$ x_{\\text{max}} =\$\$"),
            dcc_input(id="x_range2i", type="text", value="3"),
            html_label(id="y_range1", "\$\$ \\zeta(x) =\$\$"),
            dcc_input(id="y_range1i", type="text", value="-2"),
            html_label(id="y_range2", "\$\$ \\eta(x) =\$\$"),
            dcc_input(id="y_range2i", type="text", value="2"),
            html_label(id="z_range1", "\$\$ f(x, y) =\$\$"),
            dcc_input(id="z_range1i", type="text", value="-(x^2 + y^2)"),
            html_label(id="z_range2", "\$\$ g(x, y) =\$\$"),
            dcc_input(id="z_range2i", type="text", value="x^2+y^2"),
        ], style=Dict("display" => "none")
            ), html_button(id = "submit-button-state", children = "submit", n_clicks = 0),
        html_div(id="output_")
        ]),
    html_div(id="graphing", children=[
            dcc_graph(
                id = "surface_plot",
                figure = (
                    data = [
                        (x = [0],
                         y = [0],
                         z = [[0], [0]],
                        type = "surface")
                    ],
                    layout = (
                                autosize=false,
                                width=600,
                                height=600
                            )
                )
            ),
        ]
    ),
]),
html_div(id="plot_attributes",
    children=[
        html_label(id="field_density_lab", "\$\$\\text{Choose the density of the vector field.}\$\$"),
        dcc_slider(
            id = "field_density",
            min = 0,
            max = ,10
            marks = Dict([Symbol(v) => Symbol(v) for v in 0:1:10]),
            value = 5,
            step = 1,
            ),
            html_label(id="graph_accuracy_lab", "\$\$\\text{Choose the accuracy of plotting.}\$\$"),
            dcc_slider(
            id = "graph_accuracy",
            min = 10,
            max = 100,
            marks = Dict([Symbol(v) => Symbol(v) for v in 10:10:100]),
            value = 20,
            step = 10,
            )
    ]
), html_footer("\$\$\\text{MIT License}. \\ \\text{MM, ML, MK.}\$\$", style=Dict("marginTop" => "3em", "textAlign" => "center"))
])
end

callback!(app,
        Output("parametric_bounds", "style"),
        Output("fxy_bounds", "style"),
        Input("dropdown", "value")
        ) do user_choice
        if user_choice == options_[1]
            return Dict("display" => "grid"), Dict("display" => "none")
        else
            return Dict("display" => "none"), Dict("display" => "grid")
        end
end

callback!(app,
        Output("output_", "children"),
        Input("submit-button-state", "n_clicks"),
        State("dropdown", "value"),
        State("r1", "value"),
        State("r2", "value"),
        State("r3", "value"),
        State("Fx", "value"),
        State("Fy", "value"),
        State("Fz", "value"),
        State("method", "value"),
        State("integral_accuracy", "value"),
        State("u_range1i", "value"),
        State("u_range2i", "value"),
        State("v_range1i", "value"),
        State("v_range2i", "value"),
        State("x_range1i", "value"),
        State("x_range2i", "value"),
        State("y_range1i", "value"),
        State("y_range2i", "value"),
        State("z_range1i", "value"),
        State("z_range2i", "value")
        ) do n_clicks, dropdown_value, r1, r2, r3, Fx, Fy, Fz, technique, n, u_min, u_max,
            v_min, v_max, x_min, x_max, y_min, y_max, z_min, z_max
        if dropdown_value == options_[1]
            try
                u_min, u_max = parse_num.([u_min, u_max])
                p = parse_function(join(["[" * r1, r2, r3 *"]" ], ", "), :u, :v)
                q = parse_function(join(["[" * Fx, Fy, Fz * "]"], ", "), :x, :y, :z)
                ϕ = parse_function(v_min, :u)
                ψ = parse_function(v_max, :u)
                value(f, g, a, b) = round_float(Φ(f, g, (u_min, u_max), a, b; N = Int(n ÷ 2), technique = technique), 1e-3)
                value_ = @eval abs(($value($q, $p, $ϕ, $ψ)))
                return html_h5("|Φ| = $(value_).", style=Dict("textAlign" => "center"))
            catch e
                return html_h5("\$\$\\text{Wrong input. Try again.}\$\$", style=Dict("textAlign" => "center"))
            end
        else
            try
                x_min, x_max = parse_num.([x_min, x_max])
                q = parse_function(join(["[" * Fx, Fy, Fz * "]"], ", "), :x, :y, :z)
                ϕ = parse_function(y_min, :x)
                ψ = parse_function(y_max, :x)
                ρ = parse_function(z_min, :x, :y)
                η = parse_function(z_max, :x, :y)
                value(f, a, b, c, d) = round_float(Φ(f, (x_min, x_max), a, b, c, d; N = Int(n ÷ 2), technique = technique), 1e-3)
                value_ = @eval abs(($value($q, $ρ, $η, $ϕ, $ψ)))
                return html_h5("|Φ| = $(value_).", style=Dict("textAlign" => "center"))
            catch e
                return html_h5("\$\$ \\text{Wrong input. Try again.} \$\$", style=Dict("textAlign" => "center"))
            end
        end
end

callback!(app,
        Output("surface_plot", "figure"),
        Input("submit-button-state", "n_clicks"),
        State("r1", "value"),
        State("r2", "value"),
        State("r3", "value"),
        State("dropdown", "value"),
        State("graph_accuracy", "value"),
        State("u_range1i", "value"),
        State("u_range2i", "value"),
        State("v_range1i", "value"),
        State("v_range2i", "value"),
        State("x_range1i", "value"),
        State("x_range2i", "value"),
        State("y_range1i", "value"),
        State("y_range2i", "value"),
        State("z_range1i", "value"),
        State("z_range2i", "value"),
        State("Fx", "value"),
        State("Fy", "value"),
        State("Fz","value"),
        State("field_density", "value")
        ) do n_clicks, r1, r2, r3, dropdown_value, N, u_min, u_max, v_min, v_max,
            x_min, x_max, y_min, y_max, z_min, z_max, F1, F2, F3, density

            if dropdown_value == options_[1] # _______________parametric

                    try
                        u_min = parse_num(u_min)
                        u_max = parse_num(u_max)
                        v_min = parse_function(v_min, :u)
                        v_max = parse_function(v_max, :u)
                    catch e
                        u_min = -1
                        u_max = 1
                        v_min = parse_function("0", :u)
                        v_max = parse_function("2pi", :u)
                    end

                    if false
                        return Plot(surface(;
                                x = [0],
                                y = [0],
                                z = [[0], [0]]))
                    else
                        Fx = Fy = Fz = "_"
                        X(u, v) = 0
                        Y(u, v) = 0
                        Z(u, v) = 0

                        try
                            X = parse_function(r1, :u, :v)
                            Y = parse_function(r2, :u, :v)
                            Z = parse_function(r3, :u, :v)
                            Fx = parse_function(F1, :x, :y, :z)
                            Fy = parse_function(F2, :x, :y, :z)
                            Fz = parse_function(F3, :x, :y, :z)
                            @eval ($X(-2, 3), $Y(-2, 3), $Z(-2, 3))
                            @eval (isa($X(-2, 3), Array{Float64, 1}))
                            @eval (isa($Y(-2, 3), Array{Float64, 1}))
                            @eval (isa($Z(-2, 3), Array{Float64, 1}))
                            @eval (isa($Fx(-2, -2, -2), Number))
                            @eval (isa($Fy(-2, -2, -2), Number))
                            @eval (isa($Fz(-2, -2, -2), Number))
                        catch e
                            X = parse_function("0", :u, :v)
                            Y = parse_function("0", :u, :v)
                            Z = parse_function("0", :u, :v)
                            Fx = parse_function("0", :x, :y, :z)
                            Fy = parse_function("0", :x, :y, :z)
                            Fz = parse_function("0", :x, :y, :z)
                        end

                        us = LinRange(u_min, u_max, N)
                        value_(g, u::LinRange) = g.(u)
                        value_(g, u::Number) = g(u)
                        vs = @eval (LinRange(minimum($value_($v_min, $us)), maximum($value_($v_max, $us)), $N))

                        check(u, v, f) = @eval($value_($v_min, $u) <= $v && $value_($v_max, $u) >= $v ? $f($u, $v) : NaN)

                        value(g) = check.(us', vs, g)
                        x1 = @eval ($value($X))
                        y1 = @eval ($value($Y))
                        z1 = @eval ($value($Z))

                        return @eval($graph_all($x1, $y1, $z1, $Fx, $Fy, $Fz,
                                minimum(filter(!isnan, $x1)), maximum(filter(!isnan, $x1)),
                                minimum(filter(!isnan, $y1)), maximum(filter(!isnan, $y1)),
                                minimum(filter(!isnan, $z1)), maximum(filter(!isnan, $z1)), $density))
                    end
            else  # ______________________________________________f(x,y)

                Fx = Fy = Fz = "_"

                try
                    x_min = parse_num(x_min)
                    x_max = parse_num(x_max)
                    y_min = parse_function(y_min, :x)
                    y_max = parse_function(y_max, :x)
                    z_min = parse_function(z_min, :x, :y)
                    z_max = parse_function(z_max, :x, :y)
                    Fx = parse_function(F1, :x, :y, :z)
                    Fy = parse_function(F2, :x, :y, :z)
                    Fz = parse_function(F3, :x, :y, :z)
                    @eval (isa($Fx(-2, -2, -2), Number))
                    @eval (isa($Fy(-2, -2, -2), Number))
                    @eval (isa($Fz(-2, -2, -2), Number))
                catch e
                    x_min = -2
                    y_min = parse_function("-2", :x)
                    z_min = parse_function("-2", :x, :y)
                    x_max = 2
                    y_max = parse_function("2", :x)
                    z_max = parse_function("2", :x, :y)
                    Fx = parse_function("0", :x, :y, :z)
                    Fy = parse_function("0", :x, :y, :z)
                    Fz = parse_function("0", :x, :y, :z)
                end

                if false  #(x_min > x_max) | (y_min > y_max) | (z_min > z_max)
                    return Plot(surface(;
                            x = [0],
                            y = [0],
                            z = [[0], [0]]))
                else
                    xs = LinRange(x_min, x_max, N)
                    value2_(g::Function, x::LinRange{Float64}) = g.(x)
                    value2_(g::Function, x::Number) = g(x)
                    ys = @eval (LinRange(minimum($value2_($y_min, $xs)), maximum($value2_($y_max, $xs)), $N))

                    check2(x::Number, y::Number, h::Function)::Float64 = @eval($value2_($y_min, $x) <= $y && $value2_($y_max, $x) >= $y ? $h($x, $y) : NaN)

                    value2(g) = check2.(xs', ys, g)

                    zs = @eval($value2($z_max))
                    z₀ = @eval($value2($z_min))
                    temp = deepcopy(zs)
                    zs[zs .< z₀] .= NaN
                    z₀[z₀ .> temp] .= NaN

                    return @eval($graph_all($xs, $ys, $zs, $Fx, $Fy, $Fz,
                            minimum($xs), maximum($xs),
                            minimum(filter(!isnan, $ys)), maximum(filter(!isnan, $ys)),
                            minimum(filter(!isnan, $z₀)), maximum(filter(!isnan, $zs)), $density, $z₀))
                end
            end
        end

run_server(app, "0.0.0.0")


text:
# types.jl - shared types

mutable struct PolaronSite
    x::Float64
    y::Float64
    z::Float64
end


text:
using Continuables
using Documenter

makedocs(;
    modules=[Continuables],
    authors="Stephan Sahm <stephan.sahm@gmx.de> and contributors",
    repo="https://github.com/schlichtanders/Continuables.jl/blob/{commit}{path}#L{line}",
    sitename="Continuables.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://schlichtanders.github.io/Continuables.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Manual" => "manual.md",
        "Benchmark" => "benchmark.md",
        "Library" => "library.md",
    ],
)

deploydocs(;
    repo="github.com/schlichtanders/Continuables.jl",
)


text:
# Autogenerated wrapper script for HarfBuzz_jll for i686-w64-mingw32
export libharfbuzz, libharfbuzz_gobject, libharfbuzz_subset

using Cairo_jll
using Fontconfig_jll
using FreeType2_jll
using Glib_jll
using Graphite2_jll
using Libffi_jll
JLLWrappers.@generate_wrapper_header("HarfBuzz")
JLLWrappers.@declare_library_product(libharfbuzz, "libharfbuzz-0.dll")
JLLWrappers.@declare_library_product(libharfbuzz_gobject, "libharfbuzz-gobject-0.dll")
JLLWrappers.@declare_library_product(libharfbuzz_subset, "libharfbuzz-subset-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(Cairo_jll, Fontconfig_jll, FreeType2_jll, Glib_jll, Graphite2_jll, Libffi_jll)
    JLLWrappers.@init_library_product(
        libharfbuzz,
        "bin\\libharfbuzz-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_gobject,
        "bin\\libharfbuzz-gobject-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libharfbuzz_subset,
        "bin\\libharfbuzz-subset-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -1.3634297333714107 <= q <= -0.20981772159177886)

                     add_NL_constraint(m, :(log(1 + exp(-0.25086622252953283 + -0.4158912166699471*log(1 + exp(0.18862560740304835 + -0.15744636994779526*log(1 + exp(0.5575683410422836 + 0.3104874388021095*log(1 + exp(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4]))) + 0.12330243724213741*log(1 + exp(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))) + 0.6894873557418655*log(1 + exp(0.5818927234969289 + 0.4060546195732435*log(1 + exp(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4]))) + 0.22410021012329473*log(1 + exp(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))) + 0.15556574286983027*log(1 + exp(-0.6464546506721769 + -0.3052437329487665*log(1 + exp(0.5575683410422836 + 0.3104874388021095*log(1 + exp(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4]))) + 0.12330243724213741*log(1 + exp(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))) + 0.4287563119683666*log(1 + exp(0.5818927234969289 + 0.4060546195732435*log(1 + exp(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4]))) + 0.22410021012329473*log(1 + exp(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))))) + log(1 + exp(-0.9805026843034819 + -0.19393497198753096*log(1 + exp(0.18862560740304835 + -0.15744636994779526*log(1 + exp(0.5575683410422836 + 0.3104874388021095*log(1 + exp(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4]))) + 0.12330243724213741*log(1 + exp(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))) + 0.6894873557418655*log(1 + exp(0.5818927234969289 + 0.4060546195732435*log(1 + exp(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4]))) + 0.22410021012329473*log(1 + exp(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))) + -0.01642706539433636*log(1 + exp(-0.6464546506721769 + -0.3052437329487665*log(1 + exp(0.5575683410422836 + 0.3104874388021095*log(1 + exp(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4]))) + 0.12330243724213741*log(1 + exp(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))) + 0.4287563119683666*log(1 + exp(0.5818927234969289 + 0.4060546195732435*log(1 + exp(0.8764207742833086 + -0.10444718155116828*$(x[1]) + 0.8852151451027264*$(x[2]) + -0.8060184127084176*$(x[3]) + -0.617200197422199*$(x[4]))) + 0.22410021012329473*log(1 + exp(0.934307327106453 + 0.477402449923221*$(x[1]) + -0.9131747528764835*$(x[2]) + 0.07837215353953164*$(x[3]) + 0.2564992344946111*$(x[4]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
function infinite_border!(r)
    n = 1
    while isborder(r, Nord)
        if (n % 2 == 0)
            move_to_side(West, n)
        else
            move_to_side(Ost, n)
        end
        n += 1
    end
end

max_stars_repo_name: dmitrii-a-piskunov/mirea-progs

max_issues_repo_name: dmitrii-a-piskunov/mirea-progs

max_forks_repo_name: dmitrii-a-piskunov/mirea-progs


text:
using Test
using Forecast

@testset "loess" begin
    @test_throws AssertionError loess(rand(5),rand(5); d=3)

    x = loess(rand(10),rand(10); predict = rand(20))
    @test length(x) == 20

    x = loess(rand(10), rand(10); d=1, predict= rand(20))
    @test length(x) == 20
end


text:
function test_NfOrd_Elem()
  print("Testing NfOrd/Elem ...")

  Qx, x = PolynomialRing(FlintQQ, "x")

  K1, a1 = NumberField(x^3 - 2, "a")
  O1 = EquationOrder(K1)

  # creation

  b1 = @inferred O1(2*a1^0)
  @test b1.elem_in_nf == 2*a1^0
  @test parent(b1) == O1
  @test typeof(b1) == NfOrdElem{NfOrdGen}

  b2 = @inferred O1(2)
  @test parent(b2) == O1
  @test typeof(b2) == NfOrdElem{NfOrdGen}
  @test b1 == b2

  b3 = @inferred O1(fmpz(2))
  @test parent(b3) == O1
  @test typeof(b3) == NfOrdElem{NfOrdGen}
  @test b1 == b3

  b4 = @inferred O1([2, 0, 0])
  @test parent(b4) == O1
  @test typeof(b4) == NfOrdElem{NfOrdGen}
  @test b4.has_coord
  @test b1 == b4

  b5 = @inferred O1([ZZ(2), ZZ(0), ZZ(0)])
  @test parent(b5) == O1
  @test typeof(b5) == NfOrdElem{NfOrdGen}
  @test b5.has_coord
  @test b1 == b5

  b6 = @inferred O1(2*a1^0, [ZZ(2), ZZ(0), ZZ(0)])
  @test parent(b6) == O1
  @test typeof(b6) == NfOrdElem{NfOrdGen}
  @test b6.has_coord
  @test b1 == b6

  b7 = @inferred O1()
  @test parent(b6) == O1
  @test typeof(b6) == NfOrdElem{NfOrdGen}

  # deepcopy

  b = @inferred O1(a1)
  bb = @inferred deepcopy(b)
  @test b == bb

  # element in number field

  b = @inferred elem_in_nf(b1)
  @test b == K1(2)

  # element in basis

  b = @inferred elem_in_basis(b1)
  @test b == [ ZZ(2), ZZ(0), ZZ(0) ]

  b = O1(a1//2, false)
  @test_throws ErrorException elem_in_basis(b)

  # discriminant

  b = @inferred discriminant([ O1(1), O1(2*a1), O1(4*a1^2) ])
  @test b == 64 * -108

  @test_throws ErrorException discriminant( [ O1(1) ])

  # hash

  b = @inferred hash(b1)
  @test b == hash(b2)

  # zero/one

  b = @inferred one(O1)
  @test isone(b)

  b = @inferred one(b1)
  @test isone(b)

  b = @inferred zero(O1)
  @test iszero(b)

  b = @inferred zero(b1)
  @test iszero(b)

  # unary

  b = @inferred -b1
  @test b == O1(-2)

  b = -O1([2, 0, 0])
  @test b.has_coord

  # binary

  b = @inferred O1(a1) * O1(a1^2)
  @test b == O1(2)

  b = @inferred O1(a1) + O1(a1^2)
  @test b == O1(a1 + a1^2)

  c = @inferred O1([1, 0, 0]) + O1([0, 1, 0])
  @test c.has_coord

  b = @inferred O1(a1) - O1(a1^2)
  @test b == O1(a1 - a1^2)

  c = @inferred O1([1, 0, 0]) - O1([0, 1, 0])
  @test c.has_coord

  # ad hoc binary

  b = O1(2*a1)
  c = @inferred 2 * O1(a1)
  @test b == c
  c = @inferred O1(a1) * 2
  @test b == c
  c = @inferred fmpz(2) * O1(a1)
  @test b == c
  c = @inferred O1(a1) * fmpz(2)
  @test b == c

  b = O1(2 + a1)
  c = @inferred 2 + O1(a1)
  @test b == c
  c = @inferred O1(a1) + 2
  @test b == c
  c = @inferred fmpz(2) + O1(a1)
  @test b == c
  c = @inferred O1(a1) + fmpz(2)
  @test b == c

  b = O1(2 - a1)
  c = @inferred 2 - O1(a1)
  @test b == c
  c = @inferred -(O1(a1) - 2)
  @test b == c
  c = @inferred fmpz(2) - O1(a1)
  @test b == c
  c = @inferred -(O1(a1) - fmpz(2))
  @test b == c

  # ad hoc divexact

  b = O1(2*a1)
  c = @inferred divexact(b, 2)
  @test c == O1(a1)
  c = @inferred divexact(b, fmpz(2))
  @test c == O1(a1)

  # exponentiation

  b = O1(a1)
  c = @inferred b^3
  @test c == O1(2)
  c = @inferred b^fmpz(3)
  @test c == O1(2)

  # modular reduction

  b = O1([3, 2, 2])
  c = @inferred mod(b, 2)
  @test c == O1(1)
  c = @inferred mod(b, fmpz(2))
  @test c == O1(1)

  # modular exponentiation

  b = O1(2*a1)
  c = @inferred powermod(b, fmpz(3), fmpz(3))
  @test c == O1(1)
  c = @inferred powermod(b, fmpz(3), 3)
  @test c == O1(1)
  c = @inferred powermod(b, 3, fmpz(3))
  @test c == O1(1)
  c = @inferred powermod(b, 3, 3)
  @test c == O1(1)

  # representation matrix

  b = O1(1)
  c = @inferred representation_mat(b)
  @test c == one(MatrixSpace(FlintZZ, 3, 3))
  b = O1(a1)
  c = @inferred representation_mat(b)
  @test c == ZZ[0 1 0; 0 0 1; 2 0 0]
  c = @inferred representation_mat(b, K1)
  @test c == Hecke.FakeFmpqMat(ZZ[0 1 0; 0 0 1; 2 0 0], one(FlintZZ))

  # trace
  b = O1(a1)
  c = @inferred trace(b)
  @test c == 0

  # norm
  b = O1(a1)
  c = @inferred norm(b)
  @test c == 2

  # random elements
  B = 10
  b = @inferred rand(O1, -B:B)
  for i in 1:degree(O1)
    @test -B <= elem_in_basis(b)[i] && elem_in_basis(b)[i] <= B
  end

  B = BigInt(10)
  b = @inferred rand(O1, -B:B)
  B = fmpz(B)
  for i in 1:degree(O1)
    @test -B <= elem_in_basis(b)[i] && elem_in_basis(b)[i] <= B
  end

  B = 10
  b = @inferred rand(O1, B)
  for i in 1:degree(O1)
    @test -B <= elem_in_basis(b)[i] && elem_in_basis(b)[i] <= B
  end

  B = ZZ(10)
  b = @inferred rand(O1, B)
  for i in 1:degree(O1)
    @test -B <= elem_in_basis(b)[i] && elem_in_basis(b)[i] <= B
  end

  # coversion

  b = O1(a1)
  c = @inferred K1(b)
  @test c == a1

  # minkowski map

  RR = ArbField(64)

  b = O1(a1)
  c = @inferred minkowski_map(b, 1024)
  @test overlaps(c[1], root(RR(2), 3))
  @test Hecke.radiuslttwopower(c[1], -1024)
  @test overlaps(c[2],  (-1//sqrt(RR(2))*root(RR(2), 3)))
  @test Hecke.radiuslttwopower(c[2], -1024)
  @test overlaps(c[3], (sqrt(RR(3)//RR(2))*root(RR(2), 3)))
  @test Hecke.radiuslttwopower(c[3], -1024)

  # conjugates

  CC = AcbField(64)

  b = O1(a1)
  c = @inferred conjugates_arb(b, 1024)
  @test isa(c, Array{acb, 1})
  @test overlaps(c[1], CC(root(RR(2), 3)))
  @test Hecke.radiuslttwopower(c[1], -1024)
  @test overlaps(c[2], (-CC(1)//2 + onei(CC)*sqrt(RR(3))//2)*CC(root(RR(2), 3)))
  @test Hecke.radiuslttwopower(c[1], -1024)
  @test overlaps(c[3], (-CC(1)//2 - onei(CC)*sqrt(RR(3))//2)*CC(root(RR(2), 3)))
  @test Hecke.radiuslttwopower(c[1], -1024)

  c = @inferred conjugates_arb_log(b, 1024)
  @test isa(c, Array{arb, 1})
  @test overlaps(c[1], log(RR(2))//3)
  @test Hecke.radiuslttwopower(c[1], -1024)
  @test overlaps(c[2], 2*log(RR(2))//3)
  @test Hecke.radiuslttwopower(c[2], -1024)

  # test promote rules

  b = @inferred ==(O1(1), 1)

  println("PASS")
end


text:
using ImageTransformations: center

"""
    mask_circle!(::AbstractMatrix, npix; fill=0, center=center(arr))
    mask_circle!(::AbstractArray, npix; fill=0, center=center(arr))

In-place version of [`mask_circle`](@ref)
"""
function mask_circle!(arr::AbstractMatrix{T}, npix; fill=zero(T), center=center(arr)) where T
    @inbounds for idx in CartesianIndices(arr)
        d = sqrt((idx.I[1] - center[1])^2 + (idx.I[2] - center[2])^2)
        if d < npix
            arr[idx] = fill
        end
    end
    return arr
end

function mask_circle!(cube::AbstractArray{T,3}, npix; kwargs...) where T
    @inbounds for i in axes(cube, 3)
        slice = @view cube[:, :, i]
        mask_circle!(slice, npix; kwargs...)
    end
    return cube
end

"""
    mask_circle(::AbstractMatrix, npix; fill=0, center=center(arr))
    mask_circle(::AbstractArray, npix; fill=0, center=center(arr))

Mask the inner-circle of an image with radius `npix` with value `fill`. Note that the input type must be compatible with the fill value's type. If the input is a cube it will mask each frame individually.

# See Also
[`mask_circle!`](@ref)
"""
mask_circle(arr::AbstractMatrix, npix; kwargs...) = mask_circle!(copy(arr), npix; kwargs...)
mask_circle(cube::AbstractArray{T,3}, npix; kwargs...) where T = mask_circle!(copy(cube), npix; kwargs...)

"""
    mask_annulus!(::AbstractMatrix, npix_in, npix_out; fill=0, center=center(arr))

In-place version of [`mask_annulus`](@ref)
"""
function mask_annulus!(arr::AbstractMatrix{T}, npix_in, npix_out; fill = zero(T), center=center(arr)) where T
    @inbounds for idx in CartesianIndices(arr)
        d = sqrt((idx.I[1] - center[1])^2 + (idx.I[2] - center[2])^2)
        if npix_in ≤ d < npix_out
            arr[idx] = fill
        end
    end
    return arr
end

"""
    mask_annulus(::AbstractMatrix, npix_in, npix_out; fill=0, center=center(arr))

Mask an annular region of an image with inner-radius `npix_in`, outer-radius `npix_out` with value `fill`.
Note that the input type must be compatible with the fill value's type.

# See Also
[`mask_annulus!`](@ref)
"""
mask_annulus(arr::AbstractMatrix, npix_in, npix_out; kwargs...) = mask_annulus!(copy(arr), npix_in, npix_out; kwargs...)

"""
    get_annulus_segments(data, inner_radius, width, [nsegments];
                         theta_init=0, scale_factor=1, mode=:index,
                         center=center(data))

Returns indices or values in segments of a centered annulus.

### Modes
* `:mask` - returns a positive boolean mask for indexing
* `:value` - returns the data values indexed with the boolean mask
* `:apply` - returns the input data weighted by the boolean mask
"""
function get_annulus_segments(data, inner_radius, width, nsegments; theta_init=0, scale_factor=1, mode=:mask, center=center(data))
    cx, cy = center
    xs, ys = axes(data)
    dist = @. sqrt((xs - cx)^2 + (ys' - cy)^2)
    ϕ_stride = deg2rad(ceil(Int, 360 / nsegments))
    ϕ = @. atan(ys - cy, xs - cx)
    ϕ_rot = mod2pi.(ϕ)
    outer_radius = @. inner_radius + scale_factor * width

    masks = similar(data, Matrix{Bool}, nsegments)
    for i in 1:nsegments
        ϕ_start = deg2rad(theta_init) + (i - 1) * ϕ_stride
        ϕ_end = ϕ_start + ϕ_stride

        if ϕ_start < 2π && ϕ_end > 2π
            mask = @. (dist ≥ inner_radius) & (dist < outer_radius) & (ϕ_rot ≥ ϕ_start) & (ϕ_rot ≤ 2π) | (dist ≥ inner_radius) & (dist < outer_radius) & (ϕ_rot ≥ 0) & (ϕ_rot < ϕ_end - 2π)
        elseif ϕ_start ≥ 2π && ϕ_end > 2π
            mask = @. (dist ≥ inner_radius) & (dist < outer_radius) & (ϕ_rot ≥ ϕ_start - 2π) & (ϕ_rot < ϕ_end - 2π)
        else
            mask = @. (dist ≥ inner_radius) & (dist < outer_radius) & (ϕ_rot ≥ ϕ_start) & (ϕ_rot < ϕ_end)
        end
        @inbounds masks[i] = mask
    end

    return _convert_mask.(Val(mode), masks, (data,))
end

get_annulus_segments(data, inner_radius, width; kwargs...) = get_annulus_segments(data, inner_radius, width, 1; kwargs...) |> only

_convert_mask(::Val{:mask}, mask, data) = mask
_convert_mask(::Val{:value}, mask, data) = data[mask]
_convert_mask(::Val{:apply}, mask, data) = data .* mask


text:
@testset "$(exec)" for exec in [
    CompatHelper.UnableToParseSSHKey, CompatHelper.UnableToDetectCIService
]
    io = IOBuffer()
    message = "foobar"
    show(io, exec(message))

    @test contains(String(io.data), message)
end


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/subsolver_config/config.jl
# Contains subroutines used to set default configuration for select supported
# solvers along with routines needed to adjust tolerances to mirror tolerance
# adjustments in the global solve.
#############################################################################

function set_default_config_udf!(s, m::MOI.AbstractOptimizer, verbosity::Int)
    if verbosity > 0
        println("EAGO lacks a specialized configuration routine for the subsolver ($(MOI.get(m, MOI.SolverName())))")
        println("you selected. As a result, EAGO cannot set the subsolver tolerances based on the")
        println("absolute_tolerance, relative tolerance, and absolute_constraint_feas_tolerance")
        println("parameters passed to the EAGO optimizer. Consequently, need to ensure that the tolerances")
        println("set in the provided subsolver are appropriate (for instance if the absolute_tolerance = 1E-3")
        println("then the absolute tolerance for a subsolver should be < 1E-4 and any feasibility tolerances")
        println(" should be as conservative as the absolute_constraint_feas_tolerance). If you see this message")
        println("please submit an issue at https://github.com/PSORLab/EAGO.jl/issues/new/choose requesting")
        println("that a configuration routine be added for this subsolver.")
    end
    return
end

function set_default_config!(ext, d::GlobalOptimizer, m::MOI.AbstractOptimizer, local_solver)
    set_default_config_udf!(MOI.get(m, MOI.SolverName()), m, _verbosity(d))
end

function set_default_subsolver_config!(ext::DefaultExt, d::GlobalOptimizer,  m::T, local_solver::Bool) where T
    !_user_solver_config(d) && set_default_config!(ext, d, m, local_solver)
    MOI.set(m, MOI.Silent(), true)
    nothing
end

function set_default_config!(ext::DefaultExt, m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType}
    set_default_subsolver_config!(DefaultExt(), m, m._subsolvers.relaxed_optimizer, false)
    set_default_subsolver_config!(DefaultExt(), m, m._subsolvers.upper_optimizer, true)
end

"""
    set_default_config!

Configures subsolver tolerances based on tolerance parameters provided to
EAGO (provided that a specialized subsolver configuration routine has been
provided and `m.user_solver_config = false`).
"""
set_default_config!(ext::ExtensionType, m::GlobalOptimizer) = set_default_config!(DefaultExt(), m)
set_default_config!(m::GlobalOptimizer{R,S,Q}) where {R,S,Q<:ExtensionType} = set_default_config!(_ext(m), m)

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
using SauterSchwab3D: tetrahedron_face_order, tetrahedron_edge_order
using CompScienceMeshes
using SauterSchwab3D
using StaticArrays
using LinearAlgebra

const pI   = point(6,5,3)
const pII  = point(5,2,3)
const pIII = point(7,1,0)
const pIV  = point(3,4,0)

const qII  = point(3,8,4)
const qIII = point(7,5,6)
const qIV  = point(9,8,5)

const P_v = simplex(pI,pII,pIII)
const Q_v = simplex(qII,pI,qIII)

#=
s = SauterSchwab3D.singularity_detection(P_v,Q_v)
println(s)

I,J = reorder(s,cv)

K = tetrahedron_face_order(I)
L = tetrahedron_face_order(J)

println(I,J)
println(K,L)
=#

const P_e = simplex(pI,pII,pIII,pIV)
const Q_e = simplex(qIII,pI,pII,qII)

s = SauterSchwab3D.singularity_detection(P_e,Q_e)
println(s)

I,J = reorder(s)

K = tetrahedron_face_order(I)
L = tetrahedron_face_order(J)

println(I,J)
println(K,L)