text:
include("structs/dok_sparsik.jl")

using BloomFilters
using StaticArrays
using Distributions
using Profile

_powersof2 = SA[UInt64(1), UInt64(2), UInt64(4), UInt64(8), UInt64(16),
                UInt64(32), UInt64(64), UInt64(128), UInt64(256), UInt64(512),
                UInt64(1024), UInt64(2048), UInt64(4096), UInt64(8192), UInt64(16384),
                UInt64(32768), UInt64(65536), UInt64(131072), UInt64(262144), UInt64(524288),
                UInt64(1048576), UInt64(2097152), UInt64(4194304), UInt64(8388608), UInt64(16777216),
                UInt64(33554432), UInt64(67108864), UInt64(134217728), UInt64(268435456),
                UInt64(536870912), UInt64(1073741824), UInt64(2147483648), UInt64(4294967296),
                UInt64(8589934592), UInt64(17179869184), UInt64(34359738368), UInt64(68719476736),
                UInt64(137438953472), UInt64(274877906944), UInt64(549755813888), UInt64(1099511627776),
                UInt64(2199023255552), UInt64(4398046511104), UInt64(8796093022208), UInt64(17592186044416),
                UInt64(35184372088832), UInt64(70368744177664), UInt64(140737488355328), UInt64(281474976710656),
                UInt64(562949953421312), UInt64(1125899906842624), UInt64(2251799813685248), UInt64(4503599627370496),
                UInt64(9007199254740992), UInt64(18014398509481984), UInt64(36028797018963968), UInt64(72057594037927936),
                UInt64(144115188075855872), UInt64(288230376151711744), UInt64(576460752303423488), UInt64(1152921504606846976),
                UInt64(2305843009213693952), UInt64(4611686018427387904), UInt64(9223372036854775808)]

function prod_hashed(
    A::DOK_Sparsik{T}, B,
    Ahashes, Bhashes,
    hash2prod::Dict) where {T}

    # we want hash(x, y) == hash(y, x)
    hashpair(x, y) = x + y

    FF = base_ring(A)

    nnz_rows = Int[]

    rows = Dict{Int, Sparsik{T}}()
    for i in A.nnz_rows
        row_vals = Dict{Int, valtype(A)}()
        row_indices = Int[]

        for j in B.nnz_cols
            hp = hashpair(Ahashes[1][i], Bhashes[2][j])
            if haskey(hash2prod, hp)
                product = hash2prod[hp]
            else
                product = inner(A.rows[i], B.cols[j])
                hash2prod[hp] = product
            end

            if !iszero(product)
                push!(row_indices, j)
                row_vals[j] = product
            end
        end

        if !isempty(row_indices)
            push!(nnz_rows, i)
            rows[i] = Sparsik(size(A, 2), A.field, row_indices, row_vals)
        end
    end

    # we can return 0 instead of instantiating new object
    # as long as this method is not a part of interface
    return from_rows(
        size(A, 1), size(B, 2), A.field,
        nnz_rows, rows
    )
end

function do_something_hashed(matrices)
    n = size(first(matrices), 1)
    hashvector = from_dense([i in [1, 3, 5] ? 2^i : 0 for i in 1:n], QQ)
    hash2prod = Dict()

    matrixhashes = []

    for m in matrices
        row2hash, col2hash = Dict(), Dict()
        for idx in m.nnz_rows
            row2hash[idx] = inner(m.rows[idx], hashvector)
        end
        for idx in m.nnz_cols
            col2hash[idx] = inner(m.cols[idx], hashvector)
        end

        push!(matrixhashes, (row2hash, col2hash))
    end

    ans = []
    for (i, lhs) in enumerate(matrices)
        for (j, rhs) in enumerate(matrices)
            push!(
                ans,
                prod_hashed(
                    lhs,
                    rhs,
                    matrixhashes[i],
                    matrixhashes[j],
                    hash2prod)
            )
        end
    end
    ans
end

PREDICTED_ZERO = 0
ACTUAL_ZERO = 0

function fill_bloom(A::DOK_Sparsik, powersof2::Array{UInt64})
    matrixbloom = Dict{Int, UInt64}[
        Dict{Int, UInt64}(),
        Dict{Int, UInt64}()
    ]

    for idx in A.nnz_rows
        bloom = UInt64(0)
        for i in A.rows[idx].nonzero
            bloom |= powersof2[max(1, i % 64)]
        end
        matrixbloom[1][idx] = bloom
    end

    for idx in A.nnz_cols
        bloom = UInt64(0)
        for i in A.cols[idx].nonzero
            bloom |= powersof2[max(1, i % 64)]
        end
        matrixbloom[2][idx] = bloom
    end
    matrixbloom
end

function do_overlap(bv::BloomFilter, bu::BloomFilter)
    tmp = similar(bv.array)
    map!(&, bv.array, bv.array, bu.array)
    #=global PREDICTED_ZERO += Int(count(tmp) == 0)=#
    return count(tmp) != 0
end

function do_overlap(bv::UInt64, bu::UInt64)
    # global PREDICTED_ZERO += (bv & bu == 0)
    return (bv & bu) != 0
end

function prod_bloom(
    A::DOK_Sparsik{T}, B::DOK_Sparsik{T},
    bloomA::Dict{Int, UInt64}, bloomB::Dict{Int, UInt64}) where {T}

    FF = base_ring(A)
    nnz_rows = Int[]
    rows = Dict{Int, Sparsik{T}}()

    for i in A.nnz_rows
        row_vals = Dict{Int, valtype(A)}()
        row_indices = Int[]

        for j in B.nnz_cols
            #=if iszero(inner(A.rows[i], B.cols[j]))
                global ACTUAL_ZERO
                ACTUAL_ZERO += 1
            end=#
            if do_overlap(bloomA[i], bloomB[j])
                product = inner(A.rows[i], B.cols[j])
                if !iszero(product)
                    push!(row_indices, j)
                    row_vals[j] = product
                end
            end
        end

        if !isempty(row_indices)
            push!(nnz_rows, i)
            rows[i] = Sparsik(size(A, 2), A.field, row_indices, row_vals)
        end
    end

    return from_rows(
        size(A, 1), size(B, 2), A.field,
        nnz_rows, rows
    )
end

function do_something_bloom_uint(matrices)
    ans = []
    powersof2 = [UInt64(2)^i for i in 0:63]

    blooms = []
    for m in matrices
        matrixbloom = Dict{Int, UInt64}[
            Dict{Int, UInt64}(),
            Dict{Int, UInt64}()
        ]
        for idx in m.nnz_rows
            bloom = UInt64(0)
            for i in m.rows[idx].nonzero
                bloom |= powersof2[max(1, i % 64)]
            end
            matrixbloom[1][idx] = bloom
        end
        for idx in m.nnz_cols
            bloom = UInt64(0)
            for i in m.cols[idx].nonzero
                bloom |= powersof2[max(1, i % 64)]
            end
            matrixbloom[2][idx] = bloom
        end
        push!(blooms, matrixbloom)
    end

    global justbloom
    @savetime for (i, lhs) in enumerate(matrices)
        for (j, rhs) in enumerate(matrices)
            push!(ans, prod_bloom(lhs, rhs, blooms[i][1], blooms[j][2]))
        end
    end justbloom

    ans
end

function do_something_bloom_bloom(matrices)
    ans = []

    blooms = []
    @info "construct blooms"
    for m in matrices
        matrixbloom = [
            Dict{Int, BloomFilter}(),
            Dict{Int, BloomFilter}()
        ]
        for idx in m.nnz_rows
            bloom = BloomFilter(16, 0.1, 1)
            for i in m.rows[idx].nonzero
                add!(bloom, i)
            end
            matrixbloom[1][idx] = bloom
        end
        for idx in m.nnz_cols
            bloom = BloomFilter(16, 0.1, 1)
            for i in m.cols[idx].nonzero
                add!(bloom, i)
            end
            matrixbloom[2][idx] = bloom
        end
        push!(blooms, matrixbloom)
    end

    for (i, lhs) in enumerate(matrices)
        for (j, rhs) in enumerate(matrices)
            push!(ans, prod_bloom(lhs, rhs, blooms[i], blooms[j]))
            global DENS
            push!(DENS, density(last(ans)))
        end
    end

    ans
end

function hashmatrix(matrix, bloom)
    a = [UInt64(0) for _ in 1:64]
    for i in matrix.nnz_rows
        a[max(i%64, 1)] = bloom[1][i]
    end
    b = [UInt64(0) for _ in 1:64]
    for i in matrix.nnz_cols
        b[max(i%64, 1)] = bloom[2][i]
    end
    [a, b]
end

function hashmatrix_2(matrix, bloom)
    a = UInt64(0)
    for i in matrix.nnz_rows
        a |= bloom[1][i]
    end
    b = UInt64(0)
    for i in matrix.nnz_cols
        b |= bloom[2][i]
    end
    [a, b]
end

function hashmatrix_3(matrix, bloom)
    a = [UInt64(0), UInt64(0), UInt64(0), UInt64(0)]
    for i in matrix.nnz_rows
        a[max((i-1)%4+1)] |= bloom[1][i]
    end
    b = [UInt64(0), UInt64(0), UInt64(0), UInt64(0)]
    for i in matrix.nnz_cols
        b[max((i-1)%4+1)] |= bloom[2][i]
    end
    [a, b]
end

function do_overlap(bv::Array{UInt64}, bu::Array{UInt64})
    for i in length(bv)
        for j in length(bu)
            if bv[i] & bu[j] > 0
                return true
            end
        end
    end
    return false
end

function do_something_very_bloom_uint(matrices)
    ans = []
    powersof2 = [UInt64(2)^i for i in 0:63]

    blooms = [
        fill_bloom(m, powersof2)
        for m in matrices
    ]
    matrices_hashes = [
        hashmatrix_3(m, blooms[i])
        for (i, m) in enumerate(matrices)
    ]

    #println(matrices)
    #println(blooms)
    #println(matrices_hashes)

    global verybloom
    @savetime for (i, lhs) in enumerate(matrices)
        for (j, rhs) in enumerate(matrices)
            # println("overlap=$(do_overlap(matrices_hashes[i][1], matrices_hashes[j][2])) iszero=$(iszero(prod_bloom(lhs, rhs, blooms[i][1], blooms[j][2])))")
            if do_overlap(matrices_hashes[i][1], matrices_hashes[j][2])
                push!(ans, prod_bloom(lhs, rhs, blooms[i][1], blooms[j][2]))
            end
        end
    end verybloom

    ans
end

function do_something(matrices)
    ans = []

    global std
    @savetime for lhs in matrices
        for rhs in matrices
            push!(ans, lhs * rhs)

        end
    end std

    ans
end

actuall = []
predictedd = []
DENStrue = []
DENS = []

verybloom = []
justbloom = []
std = []

function UWU()

    @profile for i in 1:10
        As_std = [
            reconstruct!(from_dense(floor.(Int, abs.(randn(200, 200)) / 2), QQ))
            for _ in 1:10
        ]
        As_verysparse = [
            reconstruct!(from_dense(floor.(Int, abs.(randn(200, 200)) / 3), QQ))
            for _ in 1:10
        ]

        As = [
            As_std As_verysparse
        ]

        @info "stats" length= length(As) size=size(first(As)) density=sum(density.(As))/length(As)

        @info "normal"
        do_something(As)

        # @info "bloom bloom"
        # @time do_something_bloom_bloom(As)

        #global PREDICTED_ZERO = 0
        #global ACTUAL_ZERO = 0

        @info "bloom uint"
        do_something_bloom_uint(As)

        @info "very bloom uint"
        do_something_very_bloom_uint(As)
        #=global actuall
        global predictedd
        push!(actuall, ACTUAL_ZERO)
        push!(predictedd, PREDICTED_ZERO)=#

        println("########")
    end
end

UWU()


text:
#------------------------------------------------------------------------------
# Module
#------------------------------------------------------------------------------
type DecodeState
    ctx::Context
    global_var_num::OrderedDict{GlobalPtr, Int}
    local_var_num::OrderedDict{LLVMPtr,Int}
    local_name_counter::Int
    named_type_num::OrderedDict{TypePtr,Int}
    types_to_define::Vector{TypePtr}
end
DecodeState(ctx::Context) =
    DecodeState(ctx,
                OrderedDict{GlobalPtr,Int}(),
                OrderedDict{LLVMPtr,Int}(),
                -1,
                OrderedDict{TypePtr,Int}(),
                TypePtr[])

add_global!(st::DecodeState, val::GlobalPtr) = begin
    if haskey(st.global_var_num, val) || !isempty(FFI.get_value_name(val))
        return
    end
    n = length(st.global_var_num)
    st.global_var_num[val] = n
    return
end

get_global_name(st::DecodeState, val::GlobalPtr) = begin
    name = FFI.get_value_name(val)
    if !isempty(name)
        return Ast.Name(name)
    end
    if haskey(st.global_var_num, val)
        n = st.global_var_num[val]
    else
        n = length(st.global_var_num)
        st.global_var_num[val] = n
    end
    return Ast.UnName(n)
end

get_local_name(st::DecodeState, val) = begin
    name = FFI.get_value_name(val)
    if !isempty(name)
        return Ast.Name(name)
    end
    if haskey(st.local_var_num, val)
        n = st.local_var_num[val]
    else
        n = length(st.local_var_num)
        st.local_var_num[val] = n
        st.local_name_counter = n+1
    end
    return Ast.UnName(n)
end

get_type_name(st::DecodeState, name::String) = begin
end

save_named_type(st::DecodeState, typ::TypePtr) =
    (push!(st.types_to_define, typ); return st)

take_type_to_define(st::DecodeState) =
    (pop!(st.types_to_define); return st)

decode_llvm(st::DecodeState, tptr::TypePtr) = begin
    k = FFI.get_type_kind(tptr)
    if k == TypeKind.integer
        nbits = FFI.get_int_type_width(tptr)
        return Ast.IntType(nbits)
    elseif k == TypeKind.pointer
        etyp = decode_llvm(st, FFI.get_elem_type(tptr))
        addr = Ast.AddrSpace(FFI.get_ptr_address_space(tptr))
        return Ast.PtrType(etyp, addr)
    elseif k == TypeKind.half
        return Ast.FloatType(16, Ast.IEEE())
    elseif k == TypeKind.float
        return Ast.FloatType(32, Ast.IEEE())
    elseif k == TypeKind.double
        return Ast.FloatType(64, Ast.IEEE())
    elseif k == TypeKind.func
        return Ast.FuncType(decode_llvm(st, FFI.get_return_type(tptr)),
                            [decode_llvm(st, ty) for ty in FFI.get_param_types(tptr)],
                            FFI.is_func_var_arg(tptr))
    elseif k == TypeKind.struct
        if FFI.is_literal_struct(tptr)
            packed = FFI.is_packed_struct(tptr)
            typs = decode_llvm(st, FFI.get_struct_elem_types(tptr))
            return Ast.StructType(packed, typs)
        else
            error("unimplemented")
        end
    elseif k == TypeKind.array
        len = FFI.get_array_length(tptr)
        typ = decode_llvm(st, FFI.get_elem_type(tptr))
        return Ast.ArrayType(typ, len)
    else
        error("unimplemented type kind $k")
    end
end

decode_llvm(st::DecodeState, typs::Vector{TypePtr}) = begin
    asttyps = Array(Ast.LLVMType, length(typs))
    for i = 1:length(asttyps)
        asttyps[i] = decode_llvm(st, typs[i])
    end
    return asttyps
end

decode_llvm(st::DecodeState, cptr::ConstPtr) = begin
    ftyp = FFI.llvm_typeof(cptr)
    typ  = decode_llvm(st, ftyp)
    nops = FFI.get_num_operands(cptr)

    subclass_id = FFI.get_value_subclass_id(cptr)

    if subclass_id == ValueSubclass.undef_value
        return Ast.ConstUndef(typ)

    elseif subclass_id == ValueSubclass.global_variable
        #name = get_global_name(st, FFI.isa_global_value(cptr))
        n = st.global_var_num[FFI.isa_global_value(cptr)]
        name = Ast.UnName(n)
        return Ast.ConstGlobalRef(typ, name)

    elseif subclass_id == ValueSubclass.const_int
        words  = FFI.get_constant_int_words(cptr)
        nwords = length(words)
        # this is wrong but it gets the tests to pass
        v0 = nwords == 1 ? zero(Uint64) :
             nwords == 2 ? zero(Uint128) : zero(BigInt)
        n = foldr((b, a) -> (a << 64) | b, v0, words)
        return Ast.ConstInt(typ.nbits, n)

    elseif subclass_id == ValueSubclass.const_expr
        opcode = FFI.get_const_opcode(cptr)

        if (opcode == Opcode.add || opcode == Opcode.sub ||
            opcode == Opcode.mul || opcode == Opcode.shl)
            nsw = FFI.no_signed_wrap(cptr)
            nuw = FFI.no_unsigned_wrap(cptr)
            op1 = decode_llvm(st, FFI.get_constant_operand(cptr, 1))
            op2 = decode_llvm(st, FFI.get_constant_operand(cptr, 2))
            opcode == Opcode.add && return Ast.ConstAdd(nsw, nuw, op1, op2)
            opcode == Opcode.sub && return Ast.ConstSub(nsw, nuw, op1, op2)
            opcode == Opcode.mul && return Ast.ConstMul(nsw, nuw, op1, op2)
            opcode == Opcode.shl && return Ast.ConstShl(nsw, nuw, op1, op2)

        elseif opcode == Opcode.get_element_ptr
            inbounds = FFI.get_inbounds(cptr)
            addr = decode_llvm(st, FFI.get_constant_operand(cptr, 1))
            idxs = Array(Ast.Constant, nops-1)
            for i = 2:nops
                idxs[i-1] = decode_llvm(st, FFI.get_constant_operand(cptr, i))
            end
            return Ast.ConstGetElementPtr(inbounds, addr, idxs)

        elseif opcode == Opcode.ptrtoint
            op1 = decode_llvm(st, FFI.get_constant_operand(cptr, 1))
            return Ast.ConstPtrToInt(op1, typ)

        elseif opcode == Opcode.icmp
            pred = FFI.get_icmp_predicate(cptr)
            op1  = decode_llvm(st, FFI.get_constant_operand(cptr, 1))
            op2  = decode_llvm(st, FFI.get_constant_operand(cptr, 2))
            ipred = Ast.int_cmp_pred_map[pred]
            return Ast.ConstICmp(ipred, op1, op2)
        end
        error("unimplemented constant expr opcode $opcode")

    elseif subclass_id == ValueSubclass.const_fp
        nbits, fmt = typ.nbits, typ.fmt
        words = zeros(Uint64, div((nbits - 1), 64) + 1)
        FFI.get_const_float_words!(cptr, words)
        if nbits == 16 && typeof(fmt) === Ast.IEEE
            return Ast.ConstFloat(reinterpret(Float16, uint16(words[1])))
        elseif nbits == 32 && typeof(fmt) === Ast.IEEE
            return Ast.ConstFloat(reinterpret(Float32, uint32(words[1])))
        elseif nbits == 64 && typeof(fmt) === Ast.IEEE
            return Ast.ConstFloat(reinterpret(Float64, uint64(words[1])))
        else
            error("unimplemented constant floating point type")
        end

    elseif subclass_id == ValueSubclass.const_struct
        name = isa(typ, Ast.NamedTypeRef) ? typ.val : nothing
        p = FFI.is_packed_struct(ftyp)
        n = FFI.get_num_operands(cptr)
        vals = Array(Ast.Constant, n)
        for i = 1:n
            vals[i] = decode_llvm(st,
                FFI.get_constant_operand(cptr, i))
        end
        return Ast.ConstStruct(name, p, vals)

    elseif subclass_id == ValueSubclass.const_data_array
        @assert isa(typ, Ast.ArrayType)
        vals = Array(Ast.Constant, typ.len)
        for i = 1:typ.len
            vals[i] = decode_llvm(st,
                FFI.get_const_data_seq_elem_as_const(cptr, i))
        end
        return Ast.ConstArray(typ.typ, vals)

    elseif subclass_id == ValueSubclass.const_array
        @assert isa(typ, Ast.ArrayType)
        vals = Array(Ast.Constant, typ.len)
        for i = 1:typ.len
            vals[i] = decode_llvm(st,
                FFI.get_constant_operand(cptr, i))
        end
        return Ast.ConstArray(typ.typ, vals)

    else
        error("unimplemented subclassid : $subclass_id")
    end
end

decode_llvm(buf::MemoryBufferPtr) = begin
    start = FFI.get_buffer_start(buf)
    size  = FFI.get_buffer_size(buf)
    res = Array(Uint8, size)
    unsafe_copy!(convert(Ptr{Uint8}, res), start, size)
    return res
end
decode_llvm(st::DecodeState, buf::MemoryBufferPtr) = decode_llvm(buf)

get_metadata(iptr::InstructionPtr, n::Int) = begin
    kinds = Array(Cint, n)
    nodes = Array(MDNodePtr, n)
    nn = FFI.get_metadata(iptr, kinds, nodes, n)
    nn == 0 && return Ast.InstructionMetadata[]
    nn > n  && return get_metadata(iptr, nn)
    return zip([decode_llvm(st, kd) for kd in kinds],
               [decode_llvm(st, nd) for nd in nodes])
end
get_metadata(iptr::InstructionPtr) = get_metadata(iptr, 4)

decode_llvm(st::DecodeState, iptr::InstructionPtr) = begin
    opcode = FFI.get_instr_def_opcode(iptr)
    nops = FFI.get_num_operands(iptr)
    meta = get_metadata(iptr)
    if opcode == Opcode.ret
        op = nops == 0 ? nothing :
             Ast.ConstOperand(decode_llvm(st, FFI.get_operand(iptr, 1)))
        return Ast.Ret(op, meta)
    end
    error("unimplemented instruction opcode $opcode")
end

abstract LocalVal

immutable ForwardVal <: LocalVal
    val::Types.LLVMPtr
end

immutable DefinedVal <: LocalVal
    val::Types.LLVMPtr
end

type EncodeState
    ctx::Context
    builder::BuilderPtr
    named_types::OrderedDict{Ast.LLVMName,TypePtr}
    locals::OrderedDict{Ast.LLVMName, LocalVal}
    globals::OrderedDict{Ast.LLVMName, GlobalPtr}
    allblocks::OrderedDict{(Ast.LLVMName, Ast.LLVMName), BasicBlockPtr}
    blocks::OrderedDict{Ast.LLVMName, BasicBlockPtr}
    mdnodes::OrderedDict{Ast.MetadataNodeID, MDNodePtr}
end

EncodeState(ctx::Context) = begin
    bldr = FFI.create_builder_in_ctx(ctx)
    es = EncodeState(ctx, bldr,
                     OrderedDict{Ast.LLVMName,TypePtr}(),
                     OrderedDict{Ast.LLVMName,LocalVal}(),
                     OrderedDict{Ast.LLVMName,GlobalPtr}(),
                     OrderedDict{(Ast.LLVMName,Ast.LLVMName), BasicBlockPtr}(),
                     OrderedDict{Ast.LLVMName,BasicBlockPtr}(),
                     OrderedDict{Ast.MetadataNodeID,MDNodePtr}())
    finalizer(es, (st) -> begin
        if !isnull(st.builder)
            FFI.dispose_builder(st.builder)
            st.builder = BuilderPtr(C_NULL)
        end
    end)
    return es
end

get_named_type(st::EncodeState, name::String) = get_named_type(st, Ast.Name(name))
get_named_type(st::EncodeState, name::Ast.LLVMName) = st.named_types[name]

get_local_type(st::EncodeState, name::String) = get_local_type(st, Ast.Name(name))
get_local_type(st::EncodeState, name::Ast.LLVMName) = st.locals[name]

get_global(st::EncodeState, name::String) = get_global(st, Ast.Name(name))
get_global(st::EncodeState, name::Ast.LLVMName) = st.globals[name]

define_global!(st::EncodeState, name::String, ptr::GlobalPtr) =
    (st.globals[Ast.Name(name)] = ptr; return)

define_global!(st::EncodeState, name::Ast.LLVMName, ptr::GlobalPtr) =
    (st.globals[name] = ptr; return)

define_basic_block!(st::EncodeState,
                    fname::Ast.LLVMName,
                    bname::Ast.LLVMName,
                    ptr::BasicBlockPtr) = begin
    st.allblocks[(fname, bname)] = ptr
    st.blocks[bname] = ptr
    return
end

define_basic_block!(st::EncodeState, fname, bname, ptr) =
    define_basic_block(st, isa(fname, String) ? Ast.Name(fname) : fname::Ast.LLVMName,
                           isa(bname, String) ? Ast.Name(bname) : bname::Ast.LLVMName, ptr)

define_local!(st::EncodeState, lname::Ast.Name, ptr) = begin
    if haskey(st.locals, lname)
        def = st.locals[lanme]
        if isa(def, ForwardVal)
            FFI.replace_all_uses_with(convert(ValuePtr, def.val), ptr)
        end
        return
    end
    st.locals[lname] = DefinedVal(ptr)
    return
end

with_locals(f::Function, st::EncodeState) = begin
    ol, ob = st.locals, st.blocks
    try
        st.locals = copy(st.locals)
        st.blocks = copy(st.blocks)
        f(st)
    finally
        st.locals = ol
        st.blocks = ob
    end
end

encode_llvm(es::EncodeState, s::String) = begin
    N = length(s)
    ptr = convert(Ptr{Uint8}, Base.c_malloc(sizeof(Uint8) * N + 1))
    unsafe_copy!(ptr, convert(Ptr{Uint8}, s), N)
    unsafe_store!(ptr, zero(Uint8), N+1)
    return ptr
end

encode_llvm(es::EncodeState, typ::Ast.IntType) =
    FFI.int_type_in_ctx(es.ctx, typ.nbits)

# should all types be parametric types so we can dispatch on them more easily?
encode_llvm(st::EncodeState, typ::Ast.FloatType) = begin
    nbits, fmt = typ.nbits, typ.fmt
    if nbits == 16 && typeof(fmt) === Ast.IEEE
        return FFI.half_type_in_ctx(st.ctx)
    elseif nbits == 32 && typeof(fmt) === Ast.IEEE
        return FFI.float_type_in_ctx(st.ctx)
    elseif nbits == 64 && typeof(fmt) === Ast.IEEE
        return FFI.double_type_in_ctx(st.ctx)
    elseif nbits == 80 && typeof(fmt) === Ast.DoubleExtended
        return FFI.x86_fp80_type_in_ctx(st.ctx)
    elseif nbits == 128 && typeof(fmt) === Ast.IEEE
        return FFI.fp128_type_in_ctx(st.ctx)
    elseif nbits == 128 && typeof(fmt) === Ast.PairOfFloats
        return FFI.ppc_fp128_type_in_ctx(st.ctx)
    else
        error("unsupported floating point type $typ")
    end
end

encode_llvm(st::EncodeState, ptr::Ast.PtrType) = begin
    typ  = encode_llvm(st, ptr.typ)
    addr = encode_llvm(st, ptr.addrspace)
    return FFI.ptr_type(typ, addr)
end

encode_llvm(st::EncodeState, typs::Vector{Ast.LLVMType}) = begin
    tptrs = TypePtr[]
    for ty in typs
        push!(tptrs, encode_llvm(st, ty))
    end
    return tptrs
end

encode_llvm(st::EncodeState, cnsts::Vector{Ast.Constant}) = begin
    cptrs = ConstPtr[]
    for c in cnsts
        push!(cptrs, encode_llvm(st, c))
    end
    return cptrs
end

encode_llvm(st::EncodeState, attrs::Vector{Ast.FuncAttr}) = begin
    isempty(attrs) && return uint32(0)
    bits = uint32(attrs[1])
    for i = 2:length(attrs)
        flag = Ast.func_attr_map[attrs[i]]
        bits |= uint32(flag)
    end
    return bits
end

encode_llvm(st::EncodeState, attrs::Vector{Ast.ParamAttr}) = begin
    isempty(attrs) && return uint32(0)
    bits = uint32(attrs[1])
    for i = 2:length(attrs)
        flag = Ast.param_attr_map[attrs[i]]
        bits &= uint32(flag)
    end
    return bits
end

encode_llvm(st::EncodeState, styp::Ast.StructType) = begin
    typs = encode_llvm(st, styp.typs)
    return FFI.struct_type_in_ctx(st.ctx, typs, styp.packed)
end

encode_llvm(st::EncodeState, struct::Ast.ConstStruct) = begin
    cptrs = encode_llvm(st, struct.vals)::Vector{ConstPtr}
    if is(struct.name, nothing)
        return FFI.const_struct_in_ctx(st.ctx, cptrs, struct.packed)
    end
    # lookup named type from encoding state
    typ = st.named_types[struct.name]::TypePtr
    return FFI.const_named_struct(typ, cptrs)
end

encode_llvm(st::EncodeState, typ::Ast.ArrayType) =
    FFI.array_type(encode_llvm(st, typ.typ), typ.len)

encode_llvm(st::EncodeState, typ::Ast.FuncType) =
    FFI.func_type(encode_llvm(st, typ.rettyp),
                  encode_llvm(st, typ.argtyps),
                  typ.vaargs)

encode_llvm(st::EncodeState, carr::Ast.ConstArray) =
    FFI.const_array(encode_llvm(st, carr.typ),
                    encode_llvm(st, carr.vals))

encode_llvm(st::EncodeState, undef::Ast.ConstUndef) =
    FFI.const_undef(encode_llvm(st, undef.typ))

encode_llvm(st::EncodeState, addr::Ast.AddrSpace) = addr.val

encode_llvm(st::EncodeState, name::Ast.Name)   = name.val
encode_llvm(st::EncodeState, name::Ast.UnName) = ""

encode_llvm(st::EncodeState, ::Ast.Linkage{:External}) = 0
encode_llvm(st::EncodeState, ::Ast.DefaultVisibility)  = 0

encode_llvm(st::EncodeState, c::Ast.CConvention) = CallingCov.c

encode_llvm(st::EncodeState, val::Ast.ConstInt) = begin
    v, nbits = val.val, val.nbits
    th = encode_llvm(st, Ast.IntType(val.nbits))
    words = Uint64[(v >> 64w) & 0xffffffffffffffff for w=0:div(nbits-1, 64)]
    return FFI.const_int_arbitrary_precision(th, length(words), words)
end

encode_llvm(st::EncodeState, val::Ast.ConstFloat) = begin
    v = val.val
    nbits = 8 * sizeof(v)
    words = zeros(Uint64, div((nbits - 1), 64) + 1)
    if isa(v, Float16)
        fpsem = FloatSemantics.IEEEHalf
        words[1] = convert(Uint64, reinterpret(Uint16, v))
    elseif isa(v, Float32)
        fpsem = FloatSemantics.IEEESingle
        words[1] = convert(Uint64, reinterpret(Uint32, v))
    elseif isa(v, Float64)
        fpsem = FloatSemantics.IEEEDouble
        words[1] = reinterpret(Uint64, v)
    else
        error("unimplemented")
    end
    FFI.const_float_arbitrary_precision(st.ctx, nbits, words, fpsem)
end

encode_llvm(st::EncodeState, ref::Ast.ConstGlobalRef) =
    get_global(st, ref.name)

encode_llvm(st::EncodeState, inst::Ast.ConstPtrToInt) = begin
    val = encode_llvm(st, inst.op)
    typ = encode_llvm(st, inst.typ)
    return FFI.const_ptrtoint(val, typ)
end

encode_llvm(st::EncodeState, inst::Ast.ConstGetElementPtr) = begin
    addr = encode_llvm(st, inst.addr)
    idxs = encode_llvm(st, inst.idxs)
    return inst.inbounds ? FFI.const_inbounds_getelem_ptr(addr, idxs) :
                           FFI.const_getelem_ptr(addr, idxs)
end

encode_llvm(st::EncodeState, inst::Ast.ConstAdd) = begin
    lhs = encode_llvm(st, inst.op1)
    rhs = encode_llvm(st, inst.op2)
    inst.nsw && return FFI.const_nswadd(lhs, rhs)
    inst.nuw && return FFI.const_nuwadd(lhs, rhs)
    return FFI.const_add(lhs, rhs)
end

encode_llvm(st::EncodeState, inst::Ast.ConstICmp) = begin
    if isa(inst.pred, Ast.SGE)
        pred = IntPred.sge
    else
        error("unimplemented predicate for const icmp")
    end
    lhs = encode_llvm(st, inst.op1)
    rhs = encode_llvm(st, inst.op2)
    return FFI.const_icmp(pred, lhs, rhs)
end

# Encode Operands
encode_llvm(st::EncodeState, op::Ast.ConstOperand) =
    encode_llvm(st, op.val)

# Encode Terminators
encode_llvm(st::EncodeState, term::Ast.Ret) =
    is(term.op, nothing) ? FFI.build_ret_void(st.builder) :
                           FFI.build_ret(st.builder, encode_llvm(st, term.op))

with_sm_diagnostic(f::Function) = begin
    smd = FFI.create_sm_diagnostic()
    try
        f(smd)
    finally
        FFI.dispose_sm_diagnostic(smd)
    end
end

module_from_bitcode(ctx::Context, buf::MemoryBufferPtr) = begin
    modptr = [ModulePtr(C_NULL)]
    outmsg = Ptr{Uint8}[0]
    status = FFI.parse_llvm_bitcode(ctx, buf, modptr, outmsg)
    mod = modptr[1]
    if status != zero(LLVMBool)
        if msg[1] != C_NULL
            err = bytestring(msg[1])
            c_free(msg[1])
            throw(ErrorException(err))
        else
            error("unknown error with `parse_llvm_bitcode`")
        end
    end
    @assert !isnull(mod)
    return mod
end

# llvm takes ownership of the buffer
module_from_bitcode(ctx::Context, bytes::Union(ByteString, IOBuffer)) = begin
    buf = FFI.create_mem_buffer_with_mem_range("<bytes>", bytes.data)
    return module_from_bitcode(ctx, buf)
end

module_from_bitcode(ctx::Context, name::String, bytes::Union(ByteString, IOBuffer)) = begin
    buf = FFI.create_mem_buffer_with_mem_range(name, bytes.data)
    return module_from_bitcode(ctx, buf)
end

module_from_bitcode(ctx::Context, bytes::Vector{Uint8}) = begin
    buf = FFI.create_mem_buffer_with_mem_range("<bytes>", bytes)
    return module_from_bitcode(ctx, buf)
end

module_from_bitcode(ctx::Context, name::String, bytes::Vector{Uint8}) = begin
    buf = FFI.create_mem_buffer_with_mem_range(name, bytes)
    return module_from_bitcode(ctx, buf)
end

module_from_bitcode_file(ctx::Context, path::String) = begin
    buf = FFI.create_mem_buffer_with_contents_of_file(path)
    return module_from_bitcode(ctx, buf)
end

module_to_bitcode(ctx::Context, mod::ModulePtr) = begin
    buf = FFI.with_buff_raw_ostream() do ostream
        FFI.write_llvm_bitcode(ostream, mod)
    end
    return buf
end

module_from_assembly(ctx::Context, buf::MemoryBufferPtr) = begin
    smd = FFI.create_sm_diagnostic()
    mod = FFI.parse_llvm_assembly(ctx, buf, smd)
    try
        isnull(mod) && throw(Diagnostic(smd))
        return mod
    finally
        FFI.dispose_sm_diagnostic(smd)
    end
end

# llvm takes ownership of the buffer
module_from_assembly(ctx::Context, asm::String) = begin
    buf = FFI.create_mem_buffer_with_mem_range("<string>", asm)
    return module_from_assembly(ctx, buf)
end

module_from_assembly(ctx::Context, name::String, asm::String) = begin
    buf = FFI.create_mem_buffer_with_mem_range(name, asm)
    return module_from_assembly(ctx, buf)
end

module_from_assembly_file(ctx::Context, path::String) = begin
    buf = FFI.create_mem_buffer_with_contents_of_file(path)
    return module_from_assembly(ctx, buf)
end

module_to_assembly(ctx::Context, mod::ModulePtr) = begin
    buf = FFI.with_buff_raw_ostream() do ostream
        FFI.write_llvm_assembly(ostream, mod)
    end
    return bytestring(buf)
end

# write LLVM Assembly from a 'ModulePtr' to a file
write_assembly_file(path::String, mod::ModulePtr) = begin
    FFI.with_file_raw_ostream(path, false) do ostream
        FFI.write_llvm_assembly(ostream, mod)
    end
end

# write LLVM Bitcode from a 'ModulePtr' to a file
write_bitcode_file(path::String, mod::ModulePtr) = begin
    FFI.with_file_raw_ostream(path, false) do ostream
        FFI.write_llvm_bitcode(ostream, mod)
    end
end

get_target_triple(ptr::ModulePtr) =
    (s = FFI.get_target_triple(ptr); isempty(s) ? nothing : s)

get_datalayout(ptr::ModulePtr) =
    (s = FFI.get_datalayout(ptr); isempty(s) ? nothing : parse_datalayout(s))

module_from_ast(ctx::Context, mod::Ast.Mod) = begin
    st = EncodeState(ctx)
    mod_ptr = FFI.create_module_with_name_in_ctx(mod.name, st.ctx)
    mod.layout !== nothing && FFI.set_datalayout!(mod_ptr, mod.layout)
    mod.target !== nothing && FFI.set_target_triple!(mod_ptr, mod.target)
    # phase 1 define global definitions
    for def in mod.defs
        if isa(def, Ast.TypeDef)
            error("unimplemented")
        elseif isa(def, Ast.MetadataNodeDef)
            error("unimplemented")
        elseif isa(def, Ast.ModInlineAsm)
            error("unimplemented")
        elseif isa(def, Ast.GlobalDef)
            local g = def.val
            local gptr::GlobalValuePtr
            if isa(g, Ast.GlobalVar)
                gptr = FFI.add_global_in_addr_space!(
                                        mod_ptr,
                                        encode_llvm(st, g.typ),
                                        encode_llvm(st, g.name),
                                        encode_llvm(st, g.addrspace))
                # add global to encoding state
                define_global!(st, g.name, gptr)
            elseif isa(g, Ast.GlobalAlias)
                typ  = encode_llvm(st, g.typ)
                gptr = FFI.just_add_alias(mod_ptr,
                                          encode_llvm(st, g.typ),
                                          encode_llvm(st, g.name))
            elseif isa(g, Ast.Func)
                ftyp = encode_llvm(st,
                    Ast.FuncType(g.rettyp, [p.typ for p in g.params], g.vaargs))
                f = FFI.add_function(mod_ptr, encode_llvm(st, g.name), ftyp)
                define_global!(st, g.name, f)
            end
        end
    end
    for def in mod.defs
        if isa(def, Ast.GlobalDef)
            local g = def.val
            local gptr = get_global(st, g.name)
            if isa(g, Ast.GlobalVar)
                @assert isa(gptr, GlobalValuePtr)
                FFI.set_thread_local!(gptr, g.threadlocal)
                FFI.set_unnamed_addr!(gptr, g.unamedaddr)
                FFI.set_global_constant!(gptr, g.isconst)
                if g.init !== nothing
                    FFI.set_initializer!(gptr, encode_llvm(st, g.init))
                end
                if g.section !== nothing
                    FFI.set_section!(gptr, g.section)
                end
                FFI.set_alignment!(gptr, g.alignment)
            elseif isa(g, Ast.GlobalAlias)
                @assert isa(gptr, GlobalAliasPtr)
                FFI.set_aliasee!(gptr, encode_llvm(st, g.aliasee))
            elseif isa(g, Ast.Func)
                @assert isa(gptr, FunctionPtr)
                FFI.set_func_call_cov!(gptr, encode_llvm(st, g.callingcov))
                FFI.add_func_ret_attr!(gptr, encode_llvm(st, g.retattrs))
                FFI.add_func_attr!(gptr, encode_llvm(st, g.attrs))
                if g.section != nothing
                    FFI.set_section!(gptr, g.section)
                end
                if g.alignment != nothing
                    FFI.set_alignment!(gptr, g.alignment)
                end
                if g.gcname != nothing
                    FFI.set_gc!(gptr, g.gcname)
                end
                for blk in g.blocks
                    b = FFI.append_basicblock_in_ctx(ctx, gptr, encode_llvm(st, blk.name))
                    define_basic_block!(st, g.name, blk.name, b)
                end
                with_locals(st) do st
                    # define and add func parameters
                    pptrs = FFI.get_params(gptr)
                    for p in g.params, pptr in pptrs
                        define_local!(st, g.name, pptr)
                        FFI.set_value_name!(pptr, encode_llvm(st, p.name))
                        if !isempty(p.attrs)
                            FFI.add_attribute!(pptr, encode_llvm(st, p.attrs))
                        end
                    end
                    # build basic blocks
                    for b in g.blocks
                        bptr = st.blocks[b.name]
                        FFI.position_builder_end(st.builder, bptr)
                        # encode basic block instructions
                        for inst in b.insts
                            encode_llvm(st, inst)
                        end
                        # encode basic block terminator
                        encode_llvm(st, b.term)
                    end
                end
                # error out if any forward declared locals have not been realized
                for (k, v) in st.locals
                    isa(v, ForwardVal) && error("local $k is undefined")
                end
            end
            FFI.set_linkage!(gptr, encode_llvm(st, g.linkage))
            FFI.set_visibility!(gptr, encode_llvm(st, g.visibility))
        end
    end
    return mod_ptr
end

get_func_attrs(st::DecodeState, fnptr::FunctionPtr) = begin
    attrs = Ast.ParamAttr[]
    attr = FFI.get_func_attr(fnptr)
    bool(attr & Attribute.naked) && push!(attrs, Ast.Naked())
    bool(attr & Attribute.no_return) && push!(attrs, Ast.NoReturn())
    bool(attr & Attribute.no_unwind) && push!(attrs, Ast.NoUnwind())
    bool(attr & Attribute.read_none) && push!(attrs, Ast.ReadNone())
    bool(attr & Attribute.read_only) && push!(attrs, Ast.ReadOnly())
    bool(attr & Attribute.no_inline) && push!(attrs, Ast.NoInline())
    bool(attr & Attribute.uw_table)  && push!(attrs, Ast.UWTable())
    bool(attr & Attribute.inline_hint) && push!(attrs, Ast.InlineHint())
    bool(attr & Attribute.no_red_zone) && push!(attrs, Ast.NoRedZone())
    bool(attr & Attribute.stack_protect) && push!(attrs, Ast.StackProtect())
    bool(attr & Attribute.always_inline) && push!(attrs, Ast.AlwaysInline())
    bool(attr & Attribute.returns_twice) && push!(attrs, Ast.ReturnsTwice())
    bool(attr & Attribute.stack_alignment)   && push!(attrs, Ast.StackAlignment())
    bool(attr & Attribute.stack_protect_req) && push!(attrs, Ast.StackProtectReq())
    bool(attr & Attribute.optimize_for_size) && push!(attrs, Ast.OptimizeForSize())
    bool(attr & Attribute.no_implicit_float) && push!(attrs, Ast.NoImplicitFloat())
    return attrs
end

get_param_attrs(st::DecodeState, pptr::ParamPtr) = begin
    attrs = Ast.ParamAttr[]
    attr  = FFI.get_attribute(pptr)
    bool(attr & Attribute.s_ext)      && push!(attrs, Ast.SignExt())
    bool(attr & Attribute.in_reg)     && push!(attrs, Ast.InReg())
    bool(attr & Attribute.struct_ret) && push!(attrs, Ast.SRet())
    bool(attr & Attribute.no_alias)   && push!(attrs, Ast.NoAlias())
    bool(attr & Attribute.by_val)     && push!(attrs, Ast.ByVal())
    bool(attr & Attribute.no_capture) && push!(attrs, Ast.NoCapture())
    bool(attr & Attribute.nest)       && push!(attrs, Ast.Nest())
    return attrs
end

get_call_cov(st::DecodeState, fn::FunctionPtr) = begin
    val = FFI.get_func_call_cov(fn)
    val == CallingCov.c    && return Ast.CConvention()
    val == CallingCov.fast && return Ast.FastConvention()
    val == CallingCov.cold && return Ast.ColdConvention()
    error("unhandled calling convention enum value: $val")
end

get_linkage(st::DecodeState, fn::FunctionPtr) = begin
    val = FFI.get_visibility(fn)
    val == 0 && return Ast.Linkage{:External}()
    error("unhandled linkage enum value: $val")
end

get_visibility(st::DecodeState, fn::FunctionPtr) = begin
    val = FFI.get_visibility(fn)
    val == 0 && return Ast.DefaultVisibility()
    error("unhandled visibility enum value: $val")
end

get_params(st::DecodeState, fptr::FunctionPtr) = begin
    ps = Ast.Param[]
    for pptr in FFI.get_params(fptr)
        push!(ps, Ast.Param(decode_llvm(st, FFI.llvm_typeof(pptr)),
                            get_local_name(st, pptr),
                            get_param_attrs(st, pptr)))
    end
    return ps
end

function module_to_ast(ctx::Context, mod_ptr::ModulePtr)
    # lift c++ module to Ast.Module
    @assert ctx.handle == FFI.get_module_ctx(mod_ptr)

    st = DecodeState(ctx)

    moduleid   = FFI.get_module_id(mod_ptr)
    datalayout = get_datalayout(mod_ptr)
    triple     = get_target_triple(mod_ptr)

    # Phase 1; add toplevel global definitions
    add_toplevel(g) = add_global!(st, g)
    map(add_toplevel, FFI.list(GlobalValuePtr,
                               FFI.get_first_global(mod_ptr),
                               FFI.get_next_global))
    map(add_toplevel, FFI.list(GlobalAliasPtr,
                               FFI.get_first_alias(mod_ptr),
                               FFI.get_next_alias))
    map(add_toplevel, FFI.list(FunctionPtr,
                               FFI.get_first_func(mod_ptr),
                               FFI.get_next_func))

    # Phase 2; lift c++ definitions to Ast nodes
    local defs = Ast.Definition[]
    for g in FFI.list(GlobalValuePtr,
                      FFI.get_first_global(mod_ptr),
                      FFI.get_next_global)
        name  = get_global_name(st, g)
        ptrty = decode_llvm(st, FFI.llvm_typeof(g))
        init  = FFI.get_initializer(g)
        var = Ast.GlobalVar(name,
                            FFI.get_linkage(g),
                            FFI.get_visibility(g),
                            FFI.is_thread_local(g),
                            ptrty.addrspace,
                            FFI.has_unnamed_addr(g),
                            FFI.is_global_constant(g),
                            ptrty.typ,
                            !isnull(init) ? decode_llvm(st, init) : nothing,
                            FFI.get_section(g),
                            FFI.get_alignment(g))
        push!(defs, Ast.GlobalDef(var))
    end
    for a in FFI.list(GlobalAliasPtr,
                      FFI.get_first_alias(mod_ptr),
                      FFI.get_next_alias)
        aname = get_global_name(st, a)
        var = Ast.GlobalAlias(aname,
                              FFI.get_linkage(a),
                              FFI.get_visibility(a),
                              FFI.llvm_typeof(a),
                              decode_llvm(st, FFI.get_aliasee(a)))
        push!(defs, Ast.GlobalDef(var))
    end
    for fn in FFI.list(FunctionPtr,
                       FFI.get_first_func(mod_ptr),
                       FFI.get_next_func)
        fname = get_global_name(st, fn)
        ptrty = decode_llvm(st, FFI.llvm_typeof(fn))
        bblocks = Ast.BasicBlock[]
        for bb in FFI.list(BasicBlockPtr,
                           FFI.get_first_basicblock(fn),
                           FFI.get_next_basicblock)
            name  = get_local_name(st, bb)
            iptrs = FFI.list(InstructionPtr,
                             FFI.get_first_instruction(bb),
                             FFI.get_next_instruction)
            insts = Ast.Instruction[]
            if length(iptrs) > 2
                for i=1:length(iptrs)-1
                    push!(insts, decode_llvm(st, iptrs[i]))
                end
            end
            term  = decode_llvm(st, FFI.get_basicblock_terminator(bb))
            push!(bblocks, Ast.BasicBlock(name, insts, term))
        end
        fnty = ptrty.typ
        func = Ast.Func(get_linkage(st, fn),
                        get_visibility(st, fn),
                        get_call_cov(st, fn),
                        get_func_attrs(st, fn),
                        fnty.rettyp,
                        fname,
                        get_params(st, fn),
                        fnty.vaargs,
                        get_func_attrs(st, fn),
                        FFI.get_section(fn),
                        FFI.get_alignment(fn),
                        FFI.get_gc(fn),
                        bblocks)
        push!(defs, Ast.GlobalDef(func))
    end

    # struct definitions
    # module inline asm
    # named metadata nodes
    # metadata definitions

    return Ast.Mod(moduleid, datalayout, triple, defs)
end

# Exe Engine
remove_module(eng::ExeEnginePtr, mod::ModulePtr) = FFI.remove_module(eng, mod)

get_function(eng::ExeEnginePtr, mod::ModulePtr, fname::String) = begin
    fn = FFI.get_named_function(mod, fname)
    isnull(fn) && return nothing
    ptr = FFI.get_ptr_to_global(eng, fn)
    return ptr != C_NULL ? ptr : nothing
end

with_jit(f::Function, mod; optlevel=0) = begin
    @assert optlevel >= 0
    opt = JITCompilerOpts(0, 0, false, false)
    eng = FFI.create_mcjit_compiler_for_module(mod, opt)
    try
        f(eng)
    finally
        !isnull(eng) && FFI.dispose_execution_engine(eng)
    end
end

ext: jl


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

module Interpreter

using StaticArrays: eachindex
import Base: println, print, copy #, IteratorSize

using Raytracing, FileIO, Test
using ColorTypes:RGB
using LinearAlgebra, StaticArrays
using Printf, ProgressBars
using Documenter, DocStringExtensions, JSON

export KeywordEnum, GrammarError, InputStream
export Token, KeywordToken, IdentifierToken, StringToken
export LiteralNumberToken, SymbolToken, StopToken
export read_token, skip_whitespaces_and_comments
export Scene, InputStream, parse_scene

import Raytracing.SYM_NUM
WHITESPACE = [" ", "\t", "\n", "\r"]
OPERATIONS = ["*", "/", "+", "-"]
OPEN_BRACKETS = ["{", "[", "(", "<"]
CLOSED_BRACKETS = ["}", "]", ")", ">"]
SYMBOLS = union(OPEN_BRACKETS, CLOSED_BRACKETS, [",", "="], OPERATIONS)
LETTERS = [
     'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm',
     'n', 'o', 'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z',
     'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M',
     'N', 'O', 'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z',
     '_',
]
NUMBERS = ['0', '1', '2', '3', '4', '5', '6','7', '8', '9']

"""
   isdigit(a::String) :: Bool

Return `true` if `a` is a sigle digit, `false` otherwise.
"""
function isdigit(a::String)
     !isnothing(tryparse(Int64, a)) || (return false)
     val = parse(Int64, a)
     bool = 0≤val≤9 ? true : false
     return bool
end

"""
   isdecimal(a::String) :: Bool

Return `true` if `a` is an integer number, `false` otherwise.
"""
function isdecimal(a::String)
     !isnothing(tryparse(Int64, a)) || (return false)
     val = parse(Int64, a)
     bool = val≥0 ? true : false
     return bool
end

"""
   isalpha(a::String) :: Bool

Return `true` if `a` is a string made only of the 26 english letters (capitalized
or not) and/or the underscore symbol "_" , `false` otherwise.
"""
function isalpha(a::String)
     for ch in a
          (ch ∈ LETTERS) || (return false)
     end
     !(a=="") || (return false)
     return true
end

"""
   isalnum(a::String) :: Bool

Return `true` if `a` is a string made only of the 26 english letters (capitalized
or not), the underscore symbol "_" and the 10 basic digits, `false` otherwise.
"""
function isalnum(a::String)
     for ch in a
          (ch ∈ LETTERS || ch ∈ NUMBERS) || (return false)
     end
     !(a=="") || (return false)
     return true
end

"""
     closed_bracket(open::String) :: String

Given in input a string of an open braket ($(OPEN_BRACKETS)) return
the corresponding closed one ($(CLOSED_BRACKETS)).
"""
function close_bracket(open::String)
     (open ∈ OPEN_BRACKETS) || (throw(ArgumentError("you must insert an open braket!")))

     for i in eachindex(OPEN_BRACKETS)
          (open == OPEN_BRACKETS[i]) && (return CLOSED_BRACKETS[i])
     end

     #=
     (open == "{") || (return "}")
     (open == "[") || (return "]")
     (open == "(") || (return ")")
     (open == "<") || (return ">")
     =#

     throw(ArgumentError("you must insert an open braket!"))
end

const interpreter_dir = "interpreter"
include(joinpath(interpreter_dir, "tokens.jl"))
include(joinpath(interpreter_dir, "lexer.jl"))
include(joinpath(interpreter_dir, "parser_functions.jl"))
include(joinpath(interpreter_dir, "utilities.jl"))
include(joinpath(interpreter_dir, "parse_scene.jl"))

end # module


text:
#=
  Linear regression.

  See ~/stan/linear_regression.R
"""
#        mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
# alpha  4.03    0.00 0.13  3.79  3.95  4.03  4.12  4.28  1937    1
# beta   2.00    0.00 0.00  2.00  2.00  2.00  2.00  2.00  1975    1
# sigma  0.64    0.00 0.05  0.55  0.61  0.63  0.67  0.74  2179    1
# lp__  -4.85    0.03 1.20 -7.96 -5.42 -4.55 -3.96 -3.46  1594    1
"""

  See ~/blog/linear_regression.blog
      ~/webppl/linear_regression.wppl

  Note: I have problem with this model.
       - MH(), SMC, PG,give bad (or at least not great) rhat's
       - NUTS, HMC, etc throw errors
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function linear_regression(y)
    n = length(y)

    alpha ~ Normal(2,sqrt(2))
    beta  ~ Normal(2,sqrt(2))
    sigma ~ Gamma(2,1/2)

    # mu = Array{Float64}(undef, n)
    # mu = Vector{Float64}(undef, n)
    mu = tzeros(n)
    for i in 1:n
        mu[i] = alpha + beta*i
        y[i] ~ Normal(mu[i], sigma)
        # y[i] ~ Normal(alpha + beta*i, sigma)
    end

end

# Generated data (R):
# n = 100
# x <- 1:n
# y <- 4+2*x+rnorm(n,0,0.6)
#
# We are trying to restore these values:
#    alpha: 4.0
#    beta: 2.0
#    sigma: 0.6
#
#=
y = [5.983543,7.363850,9.001728,11.478878,14.753256,16.245218,
    18.758869,19.370087,23.417882,23.014029,25.900169,27.840894,
    29.350226,33.157245,34.046848,36.086745,38.871855,40.375736,
    42.633400,44.724727,45.845138,48.671348,50.828773,51.793457,
    54.159697,56.117874,57.545261,60.188977,62.063218,63.913408,
    66.528169,68.641637,70.790472,72.934844,73.421080,76.197684,
    78.432745,80.136338,81.058943,84.328133,85.410062,89.205925,
    89.229415,92.975393,93.531872,95.862181,97.878907,99.553297,
    102.068987,103.838061,106.568636,107.375497,110.583279,112.123316,
    113.890922,116.537840,117.444612,119.756863,123.567937,124.637554,
    125.430167,127.414399,129.799403,131.492455,134.129614,136.107248,
    138.823907,140.231367,142.460227,143.936719,146.620783,147.043767,
    150.176601,152.500756,152.473417,155.145638,158.793813,161.411058,
    163.030423,163.591843,165.696041,167.556486,169.848393,172.043780,
    173.372439,175.591097,177.409328,179.357999,181.339012,184.016542,
    186.234897,187.865354,190.316685,191.390721,193.391863,196.904094,
    197.551148,200.483832,201.944348,203.700001
    ]
=#
n = 100
x = 1:100
alpha = 4
beta = 2
sigma = 0.6
# y = alpha+beta*x+rand(Normal(0,sigma))
y = alpha.+beta .* x .+rand(Normal(0,sigma),n)
# println("y:$y")
model = linear_regression(y)

num_chains = 4

# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, MH(
#                        # :alpha => Normal(2,sqrt(2)),
#                        # :beta => Normal(2,sqrt(2)),
#                        # :sigma => Gamma(2,2)
#                        ), MCMCThreads(), 40_000, num_chains)

# chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)
chains = sample(model, SMC(1_000), MCMCThreads(), 1_000, num_chains)
# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)

# Both HMC and NUTS give the following error:
# ERROR: LoadError: TaskFailedException:
# TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing,Float64,3}
# chains = sample(model, HMC(0.1,5), MCMCThreads(), 10_000, num_chains) # Error
# chains = sample(model, NUTS(0.65), MCMCThreads(), 10_000, num_chains) # Error

display(chains)
# display(plot(chains))

# gen = generated_quantities(model, chains)
# show_var_dist_pct(gen, 40)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# ╔═╡ f9030e18-e484-11ea-24e0-3132afabdf83
begin
	using CairoMakie         # device driver for static (SVG, PDF, PNG) plots
	using Chain              # like pipes but cleaner
	using DataFrames
	using MixedModels
	using MixedModelsMakie   # plots specific to mixed-effects models using Makie
end

# ╔═╡ 1357b670-e484-11ea-2275-cba956358da4
md"# The Sleepstudy: Speed - for a change ..."

# ╔═╡ 2a6a0106-e484-11ea-2d1c-65ea0604a1e3
md"""
## Background

Belenky et al. (2003) reported effects of sleep deprivation across a 14-day study of 30-to-40-year old men and women holding commercial vehicle driving licenses. Their analyses are based on a subset of tasks and ratings from very large and comprehensive test and questionnaire battery (Balkin et al., 2000).

Initially 66 subjects were assigned to one of four time-in-bed (TIB) groups with 9 hours  (22:00-07:00) of sleep augmentation or 7 hours (24:00-07:00), 5 hours (02:00-07:00), and 3 hours (04:00-0:00) of sleep restrictions per night, respectively. The final sample comprised 56 subjects. The Psychomotor Vigilance Test (PVT) measures simple reaction time to a visual stimulus, presented approximately 10 times ⁄ minute (interstimulus interval varied from 2 to 10 s in 2-s increments) for 10 min and implemented in a thumb-operated, hand-held device (Dinges and Powell 1985).

**Design**

The study comprised 2 training days (T1, T2), one day with baseline measures (B), seven days with sleep deprivation (E1 to E7), and four recovery days (R1 to R4). T1 and T2 were devoted to training on the performance tests and familiarization with study procedures. PVT baseline testing commenced on the morning of the third day (B) and testing continued for the duration of the study (E1–E7, R1–R3; no measures were taken on R4). Bed times during T, B, and R days were 8 hours (23:00-07:00).

**Test schedule within days**

The PVT (along with the Stanford Sleepiness Scale) was administered as a battery four times per day (09:00, 12:00, 15:00, and 21:00 h); the battery included other tests not reported here (see Balkin et al. 2000). The sleep latency test was administered at 09:40 and 15:30 h for all groups. Subjects in the 3- and 5-h TIB groups performed an additional battery at 00:00 h and 02:00 h to occupy their additional time awake. The PVT and SSS were administered in this battery; however, as data from the 00:00 and 02:00 h sessions were not common to all TIB groups, these data were not included in the statistical analyses reported in the paper.

**Statistical analyses**

The authors analyzed response speed, that is (1/RT)*1000 -- completely warranted according to a Box-Cox check of the current data -- with mixed-model ANOVAs using group as between- and day as within-subject factors. The ANOVA was followed up with simple tests of the design effects implemented over days for each of the four groups.

**Current data**

The current data distributed with the _RData_ collection is attributed to the 3-hour TIB group, but the means do not agree at all with those reported for this group in Belenky et al. (2003, Figure 3) where the 3-hour TIB group is also based on only 13 (not 18) subjects. Specifically, the current data show a much smaller slow-down of response speed across E1 to E7 and do not reflect the recovery during R1 to R3. The currrent data also cover only 10 not 11 days, but it looks like only R3 is missing. The closest match of the current means was with the average of the 3-hour and 7-hour TIB groups; if only males were included, this would amount to 18 subjects. (This conjecture is based only on visual inspection of graphs.)

**References**

+ Balkin, T., et al. (2000). _Effects of sleep schedules on commercial motor vehicle driver performance._  Report MC-00–133, National Technical Information Service, U.S. Department of Transportation, Springfield, VA.

+ Belenky, G., et al. (2003). Patterns of performance degradation and restoration during sleep restriction and subsequent recovery: a sleep dose-response study. _Journal of Sleep Research_, **12**, 1-12.

+  Dinges, D. F. and Powell, J. W. (1985). Microcomputer analyses of performance on a portable, simple, visual RT task during sustained operations. _Behavior Research Methods, Instrumentation, and Computers_, **17**, 652–655.

"""

# ╔═╡ bf6eab08-e484-11ea-21b3-79b2acabd22f
md"## Setup"

# ╔═╡ 44d5f548-e4ae-11ea-13c2-ddb1bbed274e
md"""
First we attach the various packages needed, define a few helper functions, read the data, and get everything in the desired shape.
"""

# ╔═╡ e09a7f5c-1b17-4d60-958c-b84807c6638e
md"## Preprocessing"

# ╔═╡ 786e3abe-e4ae-11ea-0965-6bbd681d0f72
md"""
The `sleepstudy` data are one of the datasets available with recent versions of the `MixedModels` package. We carry out some preprocessing to have the dataframe in the desired shape:
+ Capitalize random factor `Subj`
+ Compute `speed` as an alternative dependent variable from `reaction`, warranted by a 'boxcox' check of residuals.
+ Create a `GroupedDataFrame` by levels of `Subj` (the original dataframe is available as `gdf.parent`, which we name `df`)
"""

# ╔═╡ 66a6f3bc-e485-11ea-2b47-13a48b5d7e52
gdf = @chain begin
	DataFrame(MixedModels.dataset("sleepstudy"))
	rename!(:subj => :Subj, :days => :day)
	transform!(:reaction => (x -> 1000 ./ x) => :speed)
	groupby(:Subj)
end

# ╔═╡ 61076594-0398-40c9-9190-08c820c6da97
begin
	df = gdf.parent
	describe(df)
end

# ╔═╡ 44adec0b-eec0-4b85-9161-01c11525c831
md""" ## Estimates for pooled data

In the first analysis we ignore the dependency of observations due to repeated measures from the same subjects. We pool all the data and estimate the regression of 180 speed scores on the nine days of the experiment.
"""

# ╔═╡ da75674a-aca1-4aa1-8d71-5dde799f7e90
pooledcoef = simplelinreg(df.day, df.speed)  # produces a Tuple

# ╔═╡ 6d60ce3a-2ac9-447c-8a44-d185b8562a2d
md"""
## Within-subject effects

In the second analysis we estimate coefficients for each `Subj` without regard of the information available from the complete set of data. We do not "borrow strength" to adjust for differences due to between-`Subj` variability and due to being far from the population mean.

### Within-subject simple regressions

Applying `combine` to a grouped data frame like `gdf` produces a `DataFrame` with a row for each group.
The permutation `ord` provides an ordering for the groups by increasing intercept (predicted response at day 0).
"""

# ╔═╡ 651afde1-6fe6-418e-bcd4-f95c167d7305
within = combine(gdf, [:day, :speed] => simplelinreg => :coef)

# ╔═╡ 872d5023-5d98-436d-a2c0-03492d9579b2
ord = sortperm(first.(within.coef))

# ╔═╡ 071124e6-ba1c-4d83-b5fc-97e7a45f762c
md" ### Figure of within-subject data and simple regressions"

# ╔═╡ 076b7ef7-7532-43a3-811f-78eb69cf865c
begin
	labs = values(only.(keys(gdf)))[ord]       # labels for panels
	f = clevelandaxes!(Figure(resolution=(1000,750)), labs, (2, 9))
	for (axs, sdf) in zip(f.content, gdf[ord]) # iterate over the panels and groups
		scatter!(axs, sdf.day, sdf.speed)      # add the points
		coef = simplelinreg(sdf.day, sdf.speed)
		abline!(axs, first(coef), last(coef))  # add the regression line
	end
	f
end

# ╔═╡ f3ed171c-e4ab-11ea-1eac-8b0265a6c6ea
md"## Basic LMM"

# ╔═╡ 8a902898-e485-11ea-2d50-b3c44d666f9c
m1 = fit(MixedModel, @formula(speed ~ 1+day+(1+day|Subj)), df)

# ╔═╡ 06560c54-e4ac-11ea-3b36-e57ee4bc3134
md"""
This model includes fixed effects for the intercept which estimates the average speed on the baseline day of the experiment prior to sleep deprivation, and the slowing per day of sleep deprivation. In this case about -0.11/second.

The random effects represent shifts from the typical behavior for each subject.The shift in the intercept has a standard deviation of about 0.42/s.

The within-subject correlation of the random effects for intercept and slope is small, -0.18, indicating that a simpler model with a correlation parameter (CP) forced to/ assumed to be zero may be sufficient.
"""

# ╔═╡ c01a9062-e4ac-11ea-383c-471462a1e649
md"##  No correlation parameter: zcp LMM"

# ╔═╡ 3e3b227c-e4ad-11ea-3c24-fb23d26e1171
md"""
The `zerocorr` function applied to a random-effects term estimates one parameter less than LMM `m1`-- the CP is now fixed to zero.
"""

# ╔═╡ d68c3c14-e485-11ea-0b41-6d234ab4b676
m2 = fit(MixedModel, @formula(speed ~ 1+day+zerocorr(1+day|Subj)),df)

# ╔═╡ 8e066b36-e4ad-11ea-2655-47a9a8d3c031
md"""
LMM `m2` has a slghtly  lower log-likelihood than LMM `m1` but also one fewer parameters.  The likelihood-ratio test is used to compare these nested models.
"""

# ╔═╡ ec436546-e485-11ea-19b8-a15f2f8247df
MixedModels.likelihoodratiotest(m2, m1)

# ╔═╡ d0856ae0-e4ab-11ea-0016-2f170c65768f
md"""
Alternatively, the AIC, AICc, and BIC values can be compared.  They are on a scale where "smaller is better".  All three model-fit statistics prefer the zcpLMM `m2`.
"""

# ╔═╡ 9b0956f5-e02f-4235-ba24-49e75fe71169
DataFrame(dof=dof.([m2, m1]), deviance=deviance.([m2, m1]),
          AIC = aic.([m2, m1]), AICc = aicc.([m2, m1]), BIC = bic.([m2, m1]))

# ╔═╡ a17bb040-d483-479d-8bbd-36f723f67a66
md"""
## Conditional modes of the random effects

The third set of estimates are their conditional modes. They represent a compromise between their own data and the model parameters. When distributional assumptions hold, predictions based on these estimates are more accurate than either the pooled or the within-subject estimates. Here we "borrow strength" to improve the accuracy of prediction.

"""

# ╔═╡ bc413a49-d082-40aa-8c9b-e3d64904f6cf
md" ## Three types of individual response profiles"

# ╔═╡ dab9ed8c-1deb-495b-901c-59a0c7ffcb26
md"""
Here we compare the three types on individual response profiles -- pooled, within-`Subj`, and based on conditional means of random effects.
"""

# ╔═╡ 2a424308-be20-4062-9957-887e9edf6a34
only(shrinkage(m2))

# ╔═╡ 87fd5191-9ff9-4f34-a6ed-eda19c7a6ae2
md" ## Caterpillar plots (effect profiles)"

# ╔═╡ 31b8e4c1-55d7-47b6-850a-d5356490b046
caterpillar(m2)

# ╔═╡ 168bcc03-7969-4b6a-bb2c-c1d4c18a9f4f
md" ## Shrinkage plot"

# ╔═╡ d29f6377-e8ae-40dc-b084-d31e14728edd
shrinkageplot(m2)

# ╔═╡ Cell order:
# ╟─1357b670-e484-11ea-2275-cba956358da4
# ╟─2a6a0106-e484-11ea-2d1c-65ea0604a1e3
# ╟─bf6eab08-e484-11ea-21b3-79b2acabd22f
# ╟─44d5f548-e4ae-11ea-13c2-ddb1bbed274e
# ╠═f9030e18-e484-11ea-24e0-3132afabdf83
# ╟─e09a7f5c-1b17-4d60-958c-b84807c6638e
# ╟─786e3abe-e4ae-11ea-0965-6bbd681d0f72
# ╠═66a6f3bc-e485-11ea-2b47-13a48b5d7e52
# ╠═61076594-0398-40c9-9190-08c820c6da97
# ╟─44adec0b-eec0-4b85-9161-01c11525c831
# ╠═da75674a-aca1-4aa1-8d71-5dde799f7e90
# ╟─6d60ce3a-2ac9-447c-8a44-d185b8562a2d
# ╠═651afde1-6fe6-418e-bcd4-f95c167d7305
# ╠═872d5023-5d98-436d-a2c0-03492d9579b2
# ╟─071124e6-ba1c-4d83-b5fc-97e7a45f762c
# ╠═076b7ef7-7532-43a3-811f-78eb69cf865c
# ╟─f3ed171c-e4ab-11ea-1eac-8b0265a6c6ea
# ╠═8a902898-e485-11ea-2d50-b3c44d666f9c
# ╟─06560c54-e4ac-11ea-3b36-e57ee4bc3134
# ╟─c01a9062-e4ac-11ea-383c-471462a1e649
# ╟─3e3b227c-e4ad-11ea-3c24-fb23d26e1171
# ╠═d68c3c14-e485-11ea-0b41-6d234ab4b676
# ╟─8e066b36-e4ad-11ea-2655-47a9a8d3c031
# ╠═ec436546-e485-11ea-19b8-a15f2f8247df
# ╟─d0856ae0-e4ab-11ea-0016-2f170c65768f
# ╠═9b0956f5-e02f-4235-ba24-49e75fe71169
# ╟─a17bb040-d483-479d-8bbd-36f723f67a66
# ╟─bc413a49-d082-40aa-8c9b-e3d64904f6cf
# ╟─dab9ed8c-1deb-495b-901c-59a0c7ffcb26
# ╠═2a424308-be20-4062-9957-887e9edf6a34
# ╟─87fd5191-9ff9-4f34-a6ed-eda19c7a6ae2
# ╠═31b8e4c1-55d7-47b6-850a-d5356490b046
# ╟─168bcc03-7969-4b6a-bb2c-c1d4c18a9f4f
# ╠═d29f6377-e8ae-40dc-b084-d31e14728edd


text:
# utilities for applying an FFTLog transformation

struct FFTLogPlan{T, OT, AA<:AbstractArray{Complex{T},1},
                  AAR<:AbstractArray{T,1}, PT<:Plan, IPT<:Plan}
    L::T
    N::Int
    μ::OT
    q::T
    r₀::T
    k₀r₀::T
    uₘ::AA
    r::AAR
    k::AAR
    fftplan!::PT
    ifftplan!::IPT
end

function plan_fftlog(r::AA, μ, q, k₀r₀=1.0;
                     kropt=true) where {T, AA<:AbstractArray{T,1}}
    logrmin = log(first(r))
    logrmax = log(last(r))
    r₀ = exp((logrmin + logrmax)/2)
    @assert logrmin < logrmax
    N = length(r)
    L = logrmax - logrmin
    dlnr = L / (N - 1)
    if kropt
        k₀r₀ = k₀r₀_low_ringing(N, μ, q, L, k₀r₀)
    end
    k₀ = k₀r₀ / r₀
    Nhalf = N ÷ 2
    n = range(-Nhalf,Nhalf,length=N)
    k = reverse(k₀ .* exp.(n .* L / N))

    m = fftfreq(N, N)  # get indicies that go from [-N/2] to [N/2]
    uₘ_coeff = similar(r, Complex{T})
    for i in eachindex(m)
        uₘ_coeff[i] = uₘ(m[i], μ, q, dlnr, k₀r₀, N)
    end
    uₘ_coeff[N÷2+1] = real(uₘ_coeff[N÷2+1])  # eq 19
    fftplan! = plan_fft!(uₘ_coeff)
    ifftplan! = plan_ifft!(uₘ_coeff)
    return FFTLogPlan(L, N, μ, q, r₀, k₀r₀, uₘ_coeff, r, k, fftplan!, ifftplan!)
end

U_μ(μ, x) = exp(x * log(2.) - lgamma(0.5 * (μ + 1 - x)) + lgamma(0.5 * (μ + 1 + x)))
uₘ(m, μ, q, dlnr, k₀r₀, N) = (k₀r₀)^(-2π * im * m / (dlnr * N)) * U_μ(μ, q + 2π * im * m / (dlnr * N) )

function k₀r₀_low_ringing(N, μ, q, L, k₀r₀=1.0)
    # from pyfftlog
    dlnr = L / (N-1)
    xp = (μ + 1 + q) / 2
    xm = (μ + 1 - q) / 2
    y = π * im / 2 / dlnr
    zp = lgamma(xp + y)
    zm = lgamma(xm + y)
    arg = log(2 / k₀r₀) / dlnr + imag(zp + zm) / π
    return k₀r₀ * exp((arg - round(arg))* dlnr)
end

function mul!(Y, pl::FFTLogPlan, A)
    Y .= A
    Y .*= (pl.r).^(-pl.q)
    pl.fftplan! * Y
    Y .*= pl.uₘ
    pl.ifftplan! * Y
    Y .*= (pl.r).^(pl.q)
end

function ldiv!(Y, pl::FFTLogPlan, A)
    Y .= A
    Y .*= (pl.r).^(-pl.q)
    pl.fftplan! * Y
    Y ./= pl.uₘ
    pl.ifftplan! * Y
    Y .*= (pl.r).^(pl.q)
end


text:
#!/usr/bin/env julia

using Test

@testset "FiniteVolumes.jl" begin
    # Unit tests
    include("test_mesh.jl")
    include("test_flux_functions.jl")
    include("test_scheme.jl")

    # Integration tests
    include("test_cases.jl")
    include("test_interface.jl")
end


text:
"""
    might_intersect(bb, region)

Returns true if the bounding box bb *might* intersect the region, otherwise returns
false.  To be used as an optimisation to exclude a bounding box before testing
points within the bounding box.
"""
@inline might_intersect(r1::AbstractRegion, r2::AbstractRegion) = intersects(boundingbox(r1), r2)
@inline might_intersect(bb1::BoundingBox, r2::AbstractRegion) = intersects(bb1, boundingbox(r2))
@inline might_intersect(bb::BoundingBox, region::Circle) = _might_intersect_circle(bb, region)
@inline might_intersect(bb::BoundingBox, region::Cylinder) = _might_intersect_circle(bb, region)
@inline might_intersect(bb::BoundingBox, region::Sphere) = _might_intersect_circle(bb, region)

function _might_intersect_circle(bb, circle)
    if !intersects(bb, boundingbox(circle))
        return false
    end
    # Fit a circle and find the y distance from center
    if circle.y > bb.ymax
        dy = circle.y - bb.ymax
    elseif circle.y >= bb.ymin
        dy = zero(promote_type(eltype(eltype(bb)), eltype(eltype(circle))))
    else
        dy = bb.ymin - circle.y
    end

    # Fit a circle and find the x range to query, if any
    if circle.x > bb.xmax
        dx = circle.x - bb.xmax
    elseif circle.x >= bb.xmin
        dx = zero(promote_type(eltype(eltype(bb)), eltype(eltype(circle))))
    else
        dx = bb.xmin - circle.x
    end

    h² = dx * dx + dy * dy
    r² = circle.radius * circle.radius

    return h² <= r²
end


text:
# MIT License
#
# Copyright (c) 2020 Joshua E Gentry

# 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 .Pipeline

# =========================================================================================
"""
    tess_pipeline_hocr(
        pipe::TessPipeline,
        filename::AbstractString,
        fontInfo::Bool = false
    )::Bool

Generate a HOCR file from the pipeline and save it to the specified file.  Returns `false`
if there is a problem adding the HOCR generator to the output.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | pipe     |         | The pipline to collect the HORC text from.
| R | filename |         | The file to write the HORC text to.
| O | fontInfo | `false` | Should font information be included in the output?

__Details:__

If the file exists it will be overwritten.

__Examples:__

```jldoctest; filter = r"(\\s+)"
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

tess_pipeline_hocr(pipeline, "My Book.hocr")

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

for line in readlines("My Book.hocr")[1:10]
    println(line)
end

# output

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
    "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
    <head>
        <title>My First Book</title>
        <meta http-equiv="Content-Type" content="text/html;charset=utf-8"/>
        <meta name='ocr-system' content='tesseract 4.1.1' />
        <meta name='ocr-capabilities' content='ocr_page ocr_carea ocr_par ocr_line ocrx_word ocrp_wconf'/>
    </head>
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_alto`](@ref),
          [`tess_pipeline_pdf`](@ref), [`tess_pipeline_text`](@ref)
          [`tess_pipeline_tsv`](@ref)
"""
function tess_pipeline_hocr(
            pipe::TessPipeline,
            filename::AbstractString,
            fontInfo::Bool = false
        )::Bool
    local dir         = dirname(filename)
    local tmpfile, io = mktemp(dir;cleanup=false); close(io)
    local tmpname     = basename(tmpfile)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_hocr(tmpfile, fontInfo)

    if ptr == C_NULL
        rm(tmpfile)
        return false
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local task = pipeline_start_task(
        joinpath(dir, "$(tmpname).hocr"),
        filename
    )

    pipeline_add_renderer(pipe, ptr, task)

    return true
end

# =========================================================================================
"""
    tess_pipeline_hocr(
        pipe::TessPipeline,
        fontInfo::Bool = false
    )::Union{TessOutput{String}, Nothing}

Generate an HOCR file from the pipeline and save it to a string.  Returns `nothing` if
there is a problem adding the HOCR generator to the output.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | pipe     |         | The pipline to collect the HOCR text from.
| O | fontInfo | `false` | Should font information be included in the output?

__Examples:__

```jldoctest; filter = r"(\\s+)"
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

book = tess_pipeline_hocr(pipeline)

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

for line in split(book[], "\\n")[1:10]
    println(line)
end

# output

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
    "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
    <head>
        <title>My First Book</title>
        <meta http-equiv="Content-Type" content="text/html;charset=utf-8"/>
        <meta name='ocr-system' content='tesseract 4.1.1' />
        <meta name='ocr-capabilities' content='ocr_page ocr_carea ocr_par ocr_line ocrx_word ocrp_wconf'/>
    </head>
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_alto`](@ref),
          [`tess_pipeline_pdf`](@ref), [`tess_pipeline_text`](@ref)
          [`tess_pipeline_tsv`](@ref)
"""
function tess_pipeline_hocr(
            pipe::TessPipeline,
            fontInfo::Bool = false
        )::Union{TessOutput{String}, Nothing}
    local path = mktempdir(;cleanup=false)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_hocr(
        joinpath(path, "output"),
        fontInfo
    )

    if ptr == C_NULL
        rm(path;recursive = true)
        return nothing
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local out = TessOutput(String)
    local task = pipeline_start_task(path, "output.hocr", out)

    pipeline_add_renderer(pipe, ptr, task)

    return out
end

# =========================================================================================
"""
    tess_pipeline_hocr(
        func::Function,
        pipe::TessPipeline,
        fontInfo::Bool = false
    )::Bool

Generate an HOCR file from the pipeline and pass it back to the client via a callback.
Returns `false` if there is a problem adding the HOCR generator to the output.

__Arguments:__

| T | Name     | Default | Description
|:--| :------- | :------ | :----------
| R | func     |         | The function to call with the lines of text.
| R | pipe     |         | The pipline to collect the text from.
| O | fontInfo | `false` | Should font information be included in the output?

__Details:__

The text will be passed to the caller one line at a time. The "\\n" line terminator will be
included with the text.

__Examples:__

```jldoctest; filter = r"(\\s+)"
using Tesseract

# Generate some pages to load.
write("page01.tiff", sample_tiff())
write("page02.tiff", sample_tiff())
write("page03.tiff", sample_tiff())

download_languages() # Make sure we have the data files.

instance = TessInst()
pipeline = TessPipeline(instance)

count = 0
tess_pipeline_hocr(pipeline) do line
    global count
    if count < 10
        print(line)
    end
    count += 1
end

tess_run_pipeline(pipeline, "My First Book") do add
    add(pix_read("page01.tiff"), 72)
    add(pix_read("page02.tiff"), 72)
    add(pix_read("page03.tiff"), 72)
end

# output

<?xml version="1.0" encoding="UTF-8"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN"
    "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
    <head>
        <title>My First Book</title>
        <meta http-equiv="Content-Type" content="text/html;charset=utf-8"/>
        <meta name='ocr-system' content='tesseract 4.1.1' />
        <meta name='ocr-capabilities' content='ocr_page ocr_carea ocr_par ocr_line ocrx_word ocrp_wconf'/>
    </head>
true
```

See also: [`tess_run_pipeline`](@ref), [`tess_pipeline_alto`](@ref),
          [`tess_pipeline_pdf`](@ref), [`tess_pipeline_text`](@ref)
          [`tess_pipeline_tsv`](@ref)
"""
function tess_pipeline_hocr(
            func::Function,
            pipe::TessPipeline,
            fontInfo::Bool = false
        )::Bool
    local path = mktempdir(;cleanup=false)

    # -------------------------------------------------------------------------------------
    # Create the text renderer.
    local ptr = pipeline_create_hocr(
        joinpath(path, "output"),
        fontInfo
    )

    if ptr == C_NULL
        rm(path;recursive = true)
        return false
    end

    # -------------------------------------------------------------------------------------
    # Create the task to generate the output string.
    local task = pipeline_start_task(path, "output.hocr", func, "")

    pipeline_add_renderer(pipe, ptr, task)

    return true
end


text:
#######################################
#
# elp file
#
#######################################

"""
Read elp file containing sensor locations

#### Input
* `fname`: Name or path for the sfp file
* `coordinate`: Coordinate system for electrode location
* `r`: Radius for converting spherical coords

#### Output
* `elecs`: Array of electrode objects
"""
function read_elp(fname::AbstractString; coordinate=Talairach, r::Real=90)

    @info("Reading elp file = $fname")

    # Create an empty electrode set
    elecs = Electrode[]

    # Read file and match to expected file format
    file = read(fname, String)
    regexp = r"(\S+)\s+(\S+)\s+(\S+)"
    m = matchall(regexp, file)

    # Convert label to ascii and remove '
    for idx = 1:length(m)

        local_matches = match(regexp, m[idx])

        # Extract phi and theta
        phi = float(local_matches[2])
        theta = float(local_matches[3])

        # Convert to x, y, z
        x = r .* sin.(phi*(pi/180)) .* cos.(theta*(pi/180))
        y = r .* sin.(phi*(pi/180)) .* sin.(theta*(pi/180)) - 17.5
        z = r .* cos.(phi*(pi/180))

        push!(elecs, Electrode(replace(local_matches[1], "'", "" ), coordinate(x, y, z), Dict()))
    end

    @debug("Imported $(length(elecs)) electrodes")

    return elecs
end


text:
xoff = 0
yoff = 1000

@testset "Simple Canvas" begin
    init("Simple Canvas Test", 800, 600)

    @use function setup()
        @test 800 == @width()
        @test 600 == @height()

        @test (@window) isa GtkWindow
        @test (@canvas) isa GtkCanvas
    end

    @use function update()
        w, h = @width, @height

        background("black")
        sethue("white")
        x = noise(xoff) * w
        y = noise(yoff) * h
        circle(x, y, 16, :fill)

        fontsize(18)
        text("Framerate: $(@framerate) frames/s", Point(w/2, 20), halign = :center)
        text("Framecount: $(@framecount)", Point(w/2, 40), halign = :center)

        global xoff += 0.03
        global yoff += 0.03
    end

    start(async = true)
    destroy(@window)
end

@testset "Layout & Events" begin
    init("Layout Test", 800, 600)
    @layout aside(:v, 200)

    slider = @create Slider(1:100, @margin(10), @align(:center, :fill))
    entry = @create Entry(text = "Hello world", @margin(10))
    mousepos = Ref{Point}(O)

    onevent(() -> framerate(value(slider)), :value_changed, slider)

    @use function setup()
        @test 800 == @width()
        @test 600 == @height()

        @test (@window) isa GtkWindow
        @test (@canvas) isa GtkCanvas
    end

    @use function appclosed()
        # @test true == !false
        @info "App closed"
    end

    @use function mousemove(event)
        @info "Mouse move event" pos=event.pos
        mousepos[] = event.pos
    end

    @use function mousepress(event)
        @info "Mouse press event" pos=event.pos
    end

    @use function mouserelease(event)
        @info "Mouse release event" pos=event.pos
    end

    @use function mousemotion(event)
        @info "Mouse motion event" pos=event.pos
    end

    @use function keypress(event)
        @info "Key press event" key=event.keyval
    end

    @use function update()
        w, h = @width, @height
        background("black")
        sethue("white")
        x = noise(xoff) * w
        y = noise(yoff) * h
        circle(x, y, 16, :fill)

        sethue("#f1a")
        circle(mousepos[], 16, :fill)

        global xoff += 1 / value(slider)
        global yoff += 1 / value(slider)
        noLoop()
    end

    start()
    destroy(@window)
end

@testset "Sprites" begin
    @init "Sprite test" 800 600

    sprite = loadsprite(joinpath(@__DIR__, "sprite.png"))

    @use function setup()
        @info "Sprite loaded" sprite=sprite
        @test sprite isa SpriteImage
    end

    @use function update()
        background("black")
        w, h = @width, @height
        drawsprite(sprite, Point(w/2, h/2), 400:400, centered = true)
        noLoop()
    end

    start()
    destroy(@window)
end


text:
using AbstractPlotting: stack_grouped_from_to

@testset "grouped bar: stack" begin
    x1         = [1, 1,  1,  1]
    grp_dodge1 = [2, 2,  1,  1]
    grp_stack1 = [1, 2,  1,  2]
    y1         = [2, 3, -3, -2]

    x2         = [2, 2,  2,  2]
    grp_dodge2 = [3, 4,  3,  4]
    grp_stack2 = [3, 4,  3,  4]
    y2         = [2, 3, -3, -2]

    from, to = stack_grouped_from_to(grp_stack1, y1, (; x1 = x1, grp_dodge1 = grp_dodge1))
    from1 = [0.0, 2.0,  0.0, -3.0]
    to1   = [2.0, 5.0, -3.0, -5.0]
    @test from == from1
    @test to   == to1

    from, to = stack_grouped_from_to(grp_stack2, y2, (; x2 = x2, grp_dodge2 = grp_dodge2))
    from2 = [0.0,  0.0,  0.0,  0.0]
    to2   = [2.0,  3.0, -3.0, -2.0]
    @test from == from2
    @test to   == to2

    perm = [1, 4, 2, 7, 5, 3, 8, 6]
    x = [x1; x2][perm]
    y = [y1; y2][perm]
    grp_dodge = [grp_dodge1; grp_dodge2][perm]
    grp_stack = [grp_stack1; grp_stack2][perm]

    from_test = [from1; from2][perm]
    to_test = [to1; to2][perm]

    from, to = stack_grouped_from_to(grp_stack, y, (; x = x, grp_dodge = grp_dodge))
    @test from == from_test
    @test to == to_test
end


text:
function transform(y,bxp,byp)

# This function transforms the modified variable
# to the original components of the magnetic
# perturbation according to Eq (13) of B&M (2019)
# We also impose the conditions at the pole and
# interpolate the solution to give better resolution
# in latitude

# input
#  y        - grid points in meridional coordinate
#  bxp,byp  - modified magnetic perturbation
#
# output
# ynew   - updated grid points including pole at y=1
# bx,by  - original magnetic perturbations

# apply transform back to original variables
y2 = y.*y;
yc2 = -(y2 .-1.0);
f = sqrt.(yc2);
bx = bxp .* f;
by = byp ./ f;

# append grid point at the pole y=1
 ynew = copy(y);
 push!(ynew,1.0);

# append solution at the pole
  if (m > 1)
    push!(bx,0.0);
    push!(by,0.0);
  else
    bxend = bx[end];
    byend = by[end];
    push!(bx,bxend);
    push!(by,byend);
  end

return ynew,bx,by

end

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby


text:
using LinearAlgebra: norm

"""
linearspace(start, stop, N)

Generate linear space.

Generate a linear sequence of `N` numbers between  `start` and `stop` (i. e.
sequence  of number with uniform intervals inbetween).

# Example
```
julia> linearspace(2.0, 3.0, 5)
2.0:0.25:3.0
```
"""
function linearspace(start, stop, N)
    return collect(range(Float64(start), stop = Float64(stop), length = Int64(N)))
end

"""
gradedspace(start, stop, N, strength=2)

Generate graded space.

Generate a graded sequence of `N` numbers between `start` and `stop`. This
sequence corresponds to separation of adjacent numbers that increases in
proportion corresponding to the power coefficient `strength`
from start to finish.

# Example
```
julia> gradedspace(2.0, 3.0, 5)
5-element Array{Float64,1}:
2.0
2.0625
2.25
2.5625
3.0
```
"""
function gradedspace(start, stop, N, strength=2)
    N = Int64(N)
    x = range(0.0, stop = 1.0, length = N);
    x = x.^strength
    # for i = 1:strength
    #     x = cumsum(x);
    # end
    x = x./maximum(x);
    out = Float64(start) .* (1.0 .- x) .+ Float64(stop) .* x;
end


text:
###############################################################################
#
# Atmospheric pressure and partial pressure
#
###############################################################################
GM2RT25(FT) = GRAVITY(FT) * FT(0.02896968) / RK_25(FT);

"""
    atmospheric_pressure(h::FT)

Calculate the atmospheric pressure, given
- `h` elevation in `[m]`
"""
function atmospheric_pressure_ratio(
            h::FT
) where {FT<:AbstractFloat}
    return exp( -h * GM2RT25(FT) );
end

"""
    atmospheric_pressure(h::FT)

Calculate the atmospheric pressure, given
- `h` elevation in `[m]`
"""
function atmospheric_pressure(
            h::FT
) where {FT<:AbstractFloat}
    return P_ATM(FT) * atmospheric_pressure_ratio(h);
end

"""
    ppm_to_Pa(h::FT)

Convert ppm to Pa, given
- `h` elevation in `[m]`
"""
function ppm_to_Pa(
            h::FT
) where {FT<:AbstractFloat}
    return atmospheric_pressure(h) * FT(1e-6)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
type PGConnection
    conn::Ptr{PGconn}
    closed::Bool
    busy::Bool
end

immutable PGTypeOID{tid}
end

adapt(tid) =
    let X = PGTypeOID{tid}
        method_exists(adapttype, (Type{X},)) ? adapttype(X) : String
    end

adapt(tid, vals::Vector{String}) =
    let X = PGTypeOID{tid}
        if method_exists(adapttype, (Type{X},))
            adapt(X, vals)
        else
            vals
        end
    end

function adapt{N}(X::Type{PGTypeOID{N}}, vals::Vector{String})
    T = adapttype(X)
    ovals = Vector{T}(length(vals))
    for k in eachindex(vals)
        ovals[k] = adaptval(X, vals[k])
    end
    return ovals
end

adapttype(::Type{PGTypeOID{0x00000010}}) = Bool
adapttype(::Type{PGTypeOID{0x00000012}}) = Char
adapttype(::Type{PGTypeOID{0x00000014}}) = Int64
adapttype(::Type{PGTypeOID{0x00000015}}) = Int16
adapttype(::Type{PGTypeOID{0x00000016}}) = Vector{Int16}
adapttype(::Type{PGTypeOID{0x00000017}}) = Int32
adapttype(::Type{PGTypeOID{0x0000001a}}) = UInt32
adapttype(::Type{PGTypeOID{0x0000001e}}) = Vector{UInt32}
adapttype(::Type{PGTypeOID{0x000003ed}}) = Vector{Int16}

adaptval{N}(X::Type{PGTypeOID{N}}, text::String) = parse(adapttype(X), text)

adaptval(::Type{PGTypeOID{0x00000010}}, text::String) =
    text == "t" ? true : text == "f" ? false : error()

adaptval(::Type{PGTypeOID{0x00000012}}, text::String) = text[1]

adaptval(::Type{PGTypeOID{0x00000016}}, text::String) =
    Int16[parse(Int16, word) for word in split(text)]

adaptval(::Type{PGTypeOID{0x0000001e}}, text::String) =
    UInt32[parse(UInt32, word) for word in split(text)]

adaptval(::Type{PGTypeOID{0x000003ed}}, text::String) =
    Int16[parse(Int16, word) for word in split(strip(text, ['{','}']), ',')]

function getresult(res)
    ntuples = PQntuples(res)
    nfields = PQnfields(res)
    cols = Vector{Column}(nfields)
    for j = 1:nfields
        tag  = Symbol(unsafe_string(PQfname(res, Int32(j-1))))
        tid = PQftype(res, Int32(j-1))
        nnulls = 0
        for i = 1:ntuples
            if PQgetisnull(res, Int32(i-1), Int32(j-1)) == 1
                nnulls += 1
            end
        end
        cols[j] =
            if nnulls == 0
                vals = Vector{String}(ntuples)
                for i = 1:ntuples
                    val = unsafe_string(PQgetvalue(res, Int32(i-1), Int32(j-1)))
                    vals[i] = val
                end
                OptionalColumn(OneTo(ntuples+1), adapt(tid, vals))
            elseif nnulls == ntuples
                OptionalColumn(fill(1, ntuples+1), adapt(tid, String[]))
            else
                offs = Vector{Int}(ntuples+1)
                vals = Vector{String}(ntuples-nnulls)
                offs[1] = 1
                n = 1
                for i = 1:ntuples
                    if PQgetisnull(res, Int32(i-1), Int32(j-1)) == 0
                        val = unsafe_string(PQgetvalue(res, Int32(i-1), Int32(j-1)))
                        vals[n] = val
                        n += 1
                    end
                    offs[i+1] = n
                end
                OptionalColumn(offs, adapt(tid, vals))
            end
    end
    return DataVector(Int64(ntuples), cols)
end

function connect(uri)
    pgconn = PQconnectStart(uri)
    if pgconn == C_NULL
        throw(OutOfMemoryError())
    end
    if PQstatus(pgconn) == CONNECTION_BAD
        err = unsafe_string(PQerrorMessage(pgconn))
        PQfinish(pgconn)
        error("Connection to database failed: $err")
    end
    fd = RawFD(PQsocket(pgconn))
    status = PGRES_POLLING_WRITING
    while status == PGRES_POLLING_WRITING || status == PGRES_POLLING_READING
        if status == PGRES_POLLING_WRITING
            poll_fd(fd, writable=true)
        end
        if status == PGRES_POLLING_READING
            poll_fd(fd, readable=true)
        end
        status = PQconnectPoll(pgconn)
    end
    if status != PGRES_POLLING_OK
        err = unsafe_string(PQerrorMessage(pgconn))
        PQfinish(pgconn)
        error("Connection to database failed: $err")
    end
    conn = PGConnection(pgconn, false, false)
    finalizer(conn, disconnect)
    return conn
end

function execute(conn::PGConnection, sql::String)
    if conn.closed
        error("Connection is already closed")
    end
    if conn.busy
        error("Connection is busy")
    end
    conn.busy = true
    if PQsendQuery(conn.conn, sql) != 1
        conn.busy = false
        err = unsafe_string(PQerrorMessage(conn.conn))
        error("Sending query failed: $err")
    end
    fd = RawFD(PQsocket(conn.conn))
    while PQisBusy(conn.conn) != 0
        poll_fd(fd, readable=true)
        if PQconsumeInput(conn.conn) != 1
            conn.busy = false
            err = unsafe_string(PQerrorMessage(conn.conn))
            error("Consuming input failed: $err")
        end
    end
    done = false
    out = nothing
    while !done
        res = PQgetResult(conn.conn)
        if res != C_NULL
            status = PQresultStatus(res)
            if status == PGRES_TUPLES_OK
                out = getresult(res)
            elseif status != PGRES_EMPTY_QUERY && status != PGRES_COMMAND_OK
                err = unsafe_string(PQerrorMessage(conn.conn))
                PQclear(res)
                res = PQgetResult(conn.conn)
                while res != C_NULL
                    PQclear(res)
                    res = PQgetResult(conn.conn)
                end
                conn.busy = false
                error("Reading result failed: $err")
            end
            PQclear(res)
        else
            done = true
        end
    end
    conn.busy = false
    return out
end

execute(fn::Function, conn::PGConnection, sql::String) =
    foreach(fn, execute(conn, sql))

function disconnect(conn::PGConnection)
    if !conn.closed
        PQfinish(conn.conn)
    end
    conn.closed = true
    nothing
end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto